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INTRODUCTION 


The purpose of this report is to document the main programs and 
subprograms developed and used in The Ohio State University Reports 
of the Department of Geodetic Science Nos. 156 and 157 [Papo, 1971 
and Fajemirokun, 19711, The reader is assumed to have a basic 
knowledge of the reports and terminology common to the reports 
(e. g. , coordinate system definition, etc. ) so that this documentation 
may serve as an aid for further research. 

/ Basically, this report is divided into two main sections. The 
first section deals with subroutines developed for this project. These 
subroutines constitute the mainstay of the research. The second 
section describes the main programs used to exercise these sub- 
routines to create the output given in Reports Nos. 156 and 157. 

Program/'^subprogram description is limited to formulation created 
for the reports. Other documentation mentioned in the references 
includes: 

A. Jet Propulsion Laboratory Documentation which is 
described in [O' Handley, et al. , 19691. 

B. Fortran Scientific Subroutine Documentation which 
is available in IBM publications. 

C. Ohio State University Subroutines which are general 
matrix manipulation routines (addition, subtraction, 
etc.) whose operations are easily identified by their 
names (MADD, MSUB, etc.). 

Subroutine description consists of the following elements where 
applicable: 

SUBROUTINE 
CALL STATEMENT 
SUBROUTINE PURPOSE 
INPUT PARAMETERS 
OUTPUT PARAMETERS 
COMMON AREA PARAMETERS 
PROGRAM DESCRIPTION 
SUBROUTINES REQUIRED 
REFERENCES 

Program documentation includes the above elements as well as a 
flowchart of the program logic and samples of card input. 

Extensive cross-references are used to eliminate repetition in 
descriptions of programs/subprograms using identical parameters. 

To assist in locating referenced programs an index in included in 
each major section. 
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Throughout the report floating 8 byte IBM 360 computer 
words are referred to. as "extended precision" and floating/ 
integer 4 byte computer words are referred to as "single precision 
variables/integers" respectively. 

GENERAL REFERENCES : 

A. Fajemirokun, F..A, (1971). "Application of New 
Observational Systems for Selenodetic Control, " 

The Ohio State University, Department of Geodetic 
Science, Report No. 157. 

B. O'Handley, Douglas A. et al. , (1969). "JPL Develop- 
ment Ephemeris Number 69, " Technical Report 32-1456. 

C. Papo, Haim B. (1971), "Optimal Selenodetic Control, " 
The Ohio State University, Department of Geodetic 
Science, Report No. 156. 
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DESCRIPTION OF SUBPROGRAMS 
(SUBROUTINES) 
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SOBBOUTME INDEX 


NAME 

PURPOSE 

PAGE 

APP3LTB 

Calculation of the apparent libration of the moon for an 


* 

observer on earth. 

9 

CROSVE 

Forms the vector cross product, A X B = C. 

11 

DMSRAD 

Conversion of an angle from degrees, minutes and 
seconds of arc into radians. 

13 

DVDPFl 

Administers the numerical integration of the simulated 
ephemeris of the moon through the DVDQ routine. 

15 

DVDPF2 

Slightly modified version of DVDPFl which allows usage 
of the output subroutine OUTGAJ, 

19 

EKHAHD 

Evaluation of the Eulerian angles of the moon and the 
physical librations of the moon, including their time 
rates for a given epoch. 

23 

EPPnTL 

Interpolates the 18 tabulated quantities of the simulated 
ephemeris created by the subprogram FUNEPH and 
adds the result of the interpolation to the reference 



case. 

27 

EVERAE 

Interpolation by the use of Everett’s fifth order modified 
method. 

31 

FUNEPH 

i 

Computation of the time derivatives of the "perturba- 
tions” of the state vector of the moon, the Eulerian 
angles and their time rates for the moon and the earth 
in the simulated earth-moon environment. 

35 

FUNPL5 

Computation of the second time derivatives of the phys- 
ical libration angles of the real moon, the state trans- 
ition matrix and the parameter sensitivity matri.x used 
in the numerical integration of the physical librations 
of the moon. The selenodetic positions of the earth and 
sun are obtained from the JPL DE-69 tape. 

41 

FUNPL6 

Accomplishes the same purposes as EUNPLG, but the 
simulated environment is used in lieu of the DE-69 tape. 

47 


Preceding page blank 
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NAME 


PURPOSE 


PAGE 




FUNST3 


Evaluates the time derivatives of the perturbations of 
the 18 elements of the simulated ephemeris (provided 
by FUNEPH); evaluates the time derivatives of the 
state vector of a satellite in the simulated earth-moon 
environment. 


GASTIM 


Computes the Greenwich apparent siderial time at a 
given epoch, given the nutation in obliquity and longitude 
for the epoch. 


GEULAN 

INVSPE 

LADIS 

LASOLV 
LUCA 
MAT PA 


Computes the three Eulerian angles (6, ip, <o) between 
the mean ecliptic system of 1950. 0 and the average 
terrestrial system. 

Copies a layer of a three dimensional input matrix (A) 
into a two dimensional matrix (B). 

Computes the distance between an earth observatory and 
a lunar reflector at a given epoch and the zenith distance 
of the ray at the earth observatory. 

This is the basic routine for adjusting simulated laser 
ranging. 

Creates special skew matrices for differentiation of a 
rotation matrix as outlined by Lucas, 

Forms the product A ^PA where A is an lA X JA matrix 
and P is an I A vector representing a diagonal I A X I A 
matrix and the product is dimensioned JA X JA. 


49 


55 


57 

59 


61 

65 

73 


75 


MCROSS 

MEANAN 

NUTATE 

ONEM 

OPTOBS 


Computes the cross product, R, of two input vectors, 
X and Y. 

Calculates elements of the mean orbits of the sun and 
the moon for a given epoch. 


77 

79 


Computes the nutation matrix for transformation from 
the mean to the true celestial Cartesian coordinate 
system of date. 


83 


Multiplies a matrix A by a vector B, resulting in a 

vector C. 85 

Generates a bundle of rays simulating an optical obser- 
vation from a point in space exterior to the lunar sur- 
face (e. g. , the earth or a satellite) to an arra 5 ’^ of 30 
points. 87 






NAME 


PURPOSE 





OUTADJ 

Adjusts the numerically integrated physical libration 
angles to Eckhardt's theory. 

93 

OUTEPH 

Prints routine for output of the intermediate results of 
the simulated ephemeris. 

97 

OUTGAJ 

Processes simulated earth based optical observation 
of the moon for solution of lunar coordinates and phys- 
ical parameters of the moon (initial values of the phys- 
ical libration parameters and Css, iS, Cgo)- This 
subroutine is called by the DVDPF2 routine at epochs 
where optical bundles have been observed. 

99 

PLADIS 

Computation of partial derivatives of simulated laser 
distances with respect to the considered parameters. 

107 

PMAT 

Computation of the transformation matrix necessary to 
rotate an earth fixed or lunar fixed coordinate system 
to an inertially oriented system. 

113 

POLE 

Computation of the coordinates of the true pole from 
the pole of the Conventional International Origin using 
values from the International Polar Motion Service 
(IMPS) from 1958 to 1970.5. 

115 

PRECSS 

Calculates the precession matrix P and its time deri- 
vative P. P is an orthogonal matrix that transforms a 
vector from the mean equatorial system of 1950. 0 into 
the mean equatorial system of date. 

119 

PREITR 

Prepares elements for STVITR subroutine for the 
Keplerian motion of a satellite about a primary body. 

121 

PRESTV 

Generation of series of constants and a transformation 
matrix for the calculation of the state vector of a satel- 
lite in a Keplerian orbit about a primary body. (The 
state vector is created by the use of a companion sub- 
routine STVKEP). 

123 

REDPI 

To reduce a given angle 0 to the interval 0 < 0 < 2ir. 

127 

ROTATE 

Forms a 3 X 3 rotation matrix, RNA, representing a 
rotation of angle ANG (in radians) about an axis xyz, 
designated 1, 2 and 3, I'espectively. 

129 

SEFODI 

Generation of second and fourth modified differences 
for Everett interpolation of L sets of tabulated quantities, 

131 


- 7 - 




NAME 


PURPOSE 


Transformation of a satellite state vector into instan- 
taneous Keplerian elements. 

Copies a two dimensional matrix (B) into Layer L of 
a three dimensional matrix (A), 

Computation of the rotation matrix S to rotate coordi- 
nates from the true celestial system to the average 
terrestrial system. 

Calculation of a state vector of a satellite in a Keplerian 
orbit for a given epoch and in components of a fixed 
Cartesian coordinate system. 

Calculation of a state vector in a Keplerian orbit. (This 
subroutine must be used in conjunction with subroutine 
PRESTV which generates a series of coefficients in 
common storage used in the state vector computation). 

Used with either FUNPL5 or FUNPL6. Provides the 
second time derivaUves of the physical libration para- 
meters and the ©, 9 partial derivative matrices. 

Multiplies two 3X3 matrices, A andB, to form a 
3X3 product C. 
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SUBROUTINE ; APPLIB 

CALL STATEMENT ; APPLIB (OBS, EPH, PNT, OMM, OME, AL, JOB). 

SUBROUTINE PURPOSE : Calculation of the apparent libration of the moon 
for an observer on earth. 

INPUT PARAMETERS : 

A. Vectors (extended precision) . 

1. OBS (3). Cartesian coordinates of the observer in 
the average terrestrial system {in km). 

2. EPH { 3 ). Cartesian coordinates of the moon in the 

geocentric mean ecliptic system (in km), 

3. PNT ( 3 ). Cartesian coordinates of the lunar point in the 

selenodetic system for which the position angle is 
evaluated (in km). 

B. Matrices (extended precision). 

1. OMM(3, 3). The orthogonal transformation matrix 
from the selenodetic to ecliptic system. 

2. OME (3, 3). The orthogonal transformation matrix 
from the average terrestrial system to the ecliptic 
system. 

C. Scalar (integer, single precision). 

1. JOB, An integer set to zero or non zero. If set to 

zero only the apparent librations in longitude and latitude 
are computed and the libration angles in radians are 
' returned in the first two components of the AL vector. 

If set to a non zero integer the subroutine returns (in the 
three terms of the AL vector) the total apparent libration 
in longitude, latitude and position angle. 

OUTPUT PARAMETERS : 

A. Vector (extended precision). 

1, AL( 3 ). A vector containing the apparent librations in 
longitude, latitude and (depending on the value of the 
integer JOB) the position angle (in radians). 

PROGRAM DESCRIPTION : The statements in APPLIB have .been programmed 
to,follow the expressions given in [Papo, 1971] , Section 4. 5. 
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SUBROUTINES REQUIRED : 

O. S. U. Project Library; 

1. . ONEM 

2. CROSVE 


REFERENCES : Papo, Haim B. (1971). "Optimal Selenode tic Control, '* 

Reports of the Department of Geodetic Science, No. 156, 
The Ohio State University, Columbus. 


APPIJB 


SUBRCXITINE APPLIB (CBS,EPHtPNT,OHM,CHEtAL) 

C SUBROUTINE FOR CALCULATING THE APPARENT LIBRATION OF THE MOON 

C OBS ~ CARTESIAN COORDINATES OF OBSERVER I.N AVERAGE TERESTRIAU SYSTEM 

C EPH - CARTESIAN GEOCENTRIC CCQRU3MATES OF THE HOCN IN MEAN ECLIPTIC SYSTEM 

C PNT - CARTESIAN SELENCGRAPHIC COORDINATES OF POINT ON THE HUUN 

C OBS t EPH, PNT IN KILOMETERS 

C OHM - ORIENTATION MATRIX FROM SELENOGRAPHIC TO ECLIPTIC SYSTEM 

C OHE - ORIENTATION MATRIX FROM AVERAGE TERESTRIAL TO ECLIPTIC SYSTEM 

C AL - APPARENT LIBRATION IN RADIANS 

C 11) LONGITUDE (2) LATITUDE (3) POSITION ANGLE 

IMPLICIT REAL * 8 {A-H,0-Z) 

DIMENSION OBS (3) ,EPH{3) ,OKM (3,3 ) ,OHE ( 3,3 ) , AL(3) , TU 3) ,T2{3) , T3I 3) , 
PAV(3),BV(3),CVU) ,PNTt3) ,THH(3,3} 

CALL ONEM (OME,OBS,T2) 

DO 1 1=1,3 

1 T1(I)=T2 (I)-EPH(I) 

CALL DG.MTRA {QHM,TMM,3,3 ) 

CALL QNEH ITHM,T1,T2) 

R=OSCRT{T2(l)=»’‘2+T2(2)**2+T2(3)**2) 

AL ( 1 )=DATAN2 t T2 ( 2 ), T2 ( U ) 

AL(2)=DARS1N(T2(3)/R) 

CALL ONEM (0MH,PNT,T3) 

DO 2 1=1,3 
AV(I )=DME(I,31 
BV{II=OKM(I,3) ■ 

2 CV(I>=-T1U1+T3{I) 

CALL CROSVE (CV,AV,T11 
CALL CROSVE tCV,BV,T2) 

CALL CROSVE (AV,BV,T3) 

SGN=0.00 

PA=0.00 
DO 3 1=1,3 
PA=PA+T1(I )*T2(I) 

3 SGN=SGN+CV(I)=T3( I) 

T=1.00/DSQRTUTim’>*2+Tl(2)<‘*2+Tl(3J#’»2)’»(T2(IJ**2+T2(2)#*2 +T2(3) 

P#*2>) 

AL(3) =DARCOS(PA*T) 

IF{SGN.LE-,l.O-lA)GOTO 9 
AH3)=-AH3) 

9 RETURN 
END 
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SUBROUTINE: CROSVE 


CAI,!. STATEMENT : CROSVE (A, B, C ) 
SUBROUTINE PURPOSE : Forms the vector cross product 

A X B = C 

INPUT PARAMETERS : 

A. Vectors (extended precision) 

1. A(3) 

2. B{3) 

OUTPUT PARAMETERS : 

A. Vector (extended precision) 

1. C(3) 


PROGRAM DESCRIPTION : 

A. The diagonal elements of the D matrix are initialized to be 
zero 


dij = 0 
i=j 


B. The D matrix (non diagonal terms) are then formed 
explicitly; 


d 1 2 “ ~ S3 

d2 • — a^ 

dl3 = ^2 

dj I = - ag 

das “ 

d 1 7 ~ a I 
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C, The D matrix is then multiplied by the B vector using 
the ONEM subroutine and the product is stored in the 
C vector. 


*“ — 
Cl 




bi 

Ca 




ba 

03: 


— 82 21^ 0 


K 


SUBROUTINES REQUIRED : 

O.S.U. Project Library 
1. ONEM 


REFERENCES ; None 


CROSVE 


SUBRQUTINECROSVE(AtB,C) 

C • VECTOR CROSS PRODUCT A X E = C 
IMPLICITREAL=^8 U-H,0-Z) 
DIMFNSIQNA{5) ,S(3) tC{3) ,0(3,3) 
DATA 0(1,1) ,0(2,2 ), 0(5,3) /3*0. 00/ 
0(1,2)=-A(3) 

D(2,1)=A(3) 

D(1,3)=A{2) 

D(3,1)=~A(2) 

0(2,3)--A(l) 

D(3,2)=A(1) 

Call gn£h(d,b,c) 

RETURN 

END 
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SUBROUTINE : DMSRAD 

CALL STATEMENT : DMSRAD (N, KDP, K, ANG) 

SUBROUTINE PURPOSE : Conversion of an angle from degrees , minutes 
and seconds of arc into radians. 

INPUT PARAMETERS : 

A. Integers (single precision). 

1. N contains the sign and the number of degrees, minutes 
and whole seconds. For example: 

Angle = -15° 3' 43.1305", then N = -150343 

2. KDP indicates the number of significant digits in the 
seconds fraction. For the above example KDP = 4. 

3. K is the fraction part of seconds expressed as a positive 
integer. For the above example K = 1305. 

NOTE : The subroutine is limited to N 7 ^ 0. 

OUTPUT PARAMETERS : 

A. Scalars (extended precision). 

1. ANG is the angle in radians. 

program DESCRIPTION : The input angle is converted to seconds of arc and 
converted to radians. 


SUBROUTINES REQUIRED : None 


REFERENCES : None 
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DMSRAD 


C 

C 

c 


SU3RQUTINE OMSRAD (N , KGP ,K , ANG) 

CONVERSION FROM OEG/HIN/SEC TO RADIANS 

IMPORTANT : N StIOULO DE DIFFERENT FRDH ZERO 

N - SIGN , DEGREES MINUTES AND SECONDS (INTEGERS) 

K - FRACTION OF A SECOND WITH KOP DIGITS 
IMPLICIT REAL *8 CA-H,0-Z) 

IF(N.LT.O)K=-R 

NOM=N/lCO 

ND=NDM/100 

ANG— tOPLOAT (ND^3600+i‘4H*60+N— W0»M^1C0)+DFL0AT (K)/l 10^D0**KDP) )>-48^8 

P1368U09536D-05 

RETURN 

ENa 
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SUBROUTINE : DVDPFl 


CALL STATEMENT : D VD PF 1 (Parameters as in DVDQ, RAT, RAT2, 

IHL, PAR, FUN, OUT, JOB) 

SUBROUTINE PURPOSE ; To administer the numerical integration of the simulated 
ephemeris through the use of the DVDQ numerical integrator and customer 
supplied subroutines. 

INPUT PARAMETERS : 

A. Parameters input to the DVDQ routine are described in 
fKrogh, 1969], 

B. - Scalars (extended precision). 

1. RATI is relative accuracy required in integration for the 
ephemeris of the moon (e. g. the geocentric state vector), 

2. RAT2 is the relative accuracy required in integration of 
theTEulerian angles of the earth-moon system (in radians). 

C- Integers (single precision). 

1. IHL is a dummy variable. 

2. JOB is a control integer for use by the subroutine. If 
set to zero, the subroutine will perform integration only. 

If set to 1, the subroutine will also allow the accumulation 
of normal equations and a constant vector for a least 
squares adjustment. 

D. Vector (extended precision). 

1. PAR is a vector of parameters input from the main program 
and used for generating the simulated ephemeris and descril> 
in the documentation of the main program. 

E. External subroutine. 

1. FUN is the name of the subroutine which generates time 
derivatives for a given set of functions (FUNEPH for 
example). 

OUTPUT PARAMETERS : 

A. External subroutine. 

1. OUT is the name of the output subroutine for passing control 
over to DVDPFl at selected input intervals (the parameter 
BELT). 

B. Common area DVOUT is used only if JOB = 1, and contains the 
following; 

1. Scalar (extended precision). 
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a, ETOB is the next epoch at which an obser- 
vation has been made and program control 
is to be passed to OUT. 

2. Integers (single precision). 

a. NUT is a flag which is set to zero if the 
integration is proceeding normally and to one 
if errors are detected in FUN or OUT. If set 
to one, the integration is terminated. 

b. ICW is the serial number of the next obseirva- 
tion. 

PROGRAM DESCRIPTION : The program is the primary ''integrator” of the parameters 
in the simulated environment fPapo, 19711. The basic integrator is 
DVDQ and the calculation of time derivatives is provided by FUNEPH. 
Description of the DVDQ is contained in fKrogh, 1969] and the 
equations assumed for integration in [ Papo,. 1971] Chapter 4. 


SUBROUTINES REQUIRED : 

O. S.U. Project Library: 

1. DVDQ (JPL subroutine) 

2. DVDQl (JPL subroutine altered for input after initial call 
of DVDQ) 

3. External subroutines given in input-output section. 


REFERENCES : 

A. Krogh, F. T. (1969). "VODQ/SVDQ/DVDQ— Variable 
Order Integrators for the Numerical Solution of Ordinary 
Differential Equations, " Section 3. 14, JPL, Pasadena, 
California. 

B. Papo, Haim B. (1971). "Optimal Selenodetic Control, ” 
The Ohio State University, Department of Geodetic Science , 
Report No. 156 . 
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non 


DVDPFl 


SUBROUTINE OVOPFKNQ , Y , YP , KD . £ P . RATI ,RATZ ,HHlf^HMAX , DELT 

P »KST,lHLtPAR,FUNfOUT?KUTYN,CTf JOS] 

SUBROUTINE FOR INTEGRATING OIFERENI lAL EQUATIONS BY THE DVUQ 
NUMERICAL INTEGRATOR . THE NUMBER UF INTEGRANDS SHUULO BE AT 
LEAST SIX , 

IMPLICIT REALMS (A-H,0-2) 

REAL*4 EP(NQ ) »HMIN,HMAX, EHX,RAT1 ,RAT2 

HXSTP=500000 
K0C1}=1 
ET=PARU ) 

ETFN=PAR(2) 

H=PAR(3) 

NUT=0 

1CW=0 

HRITE(6,8 9UY,ET,ETFN,H,DELT 

HRITE(6,892 iHMlNtHHAX »RAT1 , RaT2 ,NQ,KO 1 1 ) »HXSTP 

KLU=-1 

GOTO 881 

879 CALL DVOC tNQ,ET, Y»YPtKDtEP»lFLiH»HHIN»HMAX,DELT,ETFN,HXSTP, 
PKST.K MX,EMX.KQ,YN,OT] 

GOTO 884 

863 CALL DVDQl (NQ tE T» Y, yp ,KO , EP t IFL tH tUMlN fHMAX»DELT»fcTFNtMXSTP, 
PKST.K MX.EHX.KO.YN.OT) 

884 GOTO (881, 882, S85.685»8e5»686. 887.888) , IFL 

885 CALL FUN (ET,Y,YP) 

CALL 0UT(ET,Y,YP,1HL,NQ,PAR) 

IFtJCB.NE.l) goto 883 
DELT=ETOe-ET 
■ ICH=ICW+1 

IFCNUT.EO.IJ GOTO 838 
GOTO 883 

866 EP(KHX)=32.*EHX*EP(KHX) 

IF (EPIKHXJ.GT.O. ) EP (KHX)=-EP(KMX) 

GOTO 083 
881 DO 871 1=1,6 

EP(1)=Y(I)*RAT1 

1F(A8S(EP(I) ).LT.5.E-07) EP (I)=-5-E-07 

871 IF(EP(I).GT.0.) EPm=-EP(I) 

IF(N0-6 ) 874,874,872 

872 DO 873 1=7,12 
EP(I )=Y(I )*RAT2 

EP( 1+6)=Y( I+6)’»RAT2*10000. 

IF(A8S{EP(I)).LT.1.E~07) cP ( 1 J=-1 .E-07 
IF(ABStEP{I-*-6) ).LT.l.E-09) EP ( I-t-6 ) =-l .E-C9 

873 IF( EP(I ) .GT.O. » EP(I)=-EP(1> 

874 CONTINUE 
KLU=KLU+1 
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882 CALL FUN (ET.YtYP) 

IF(KLU) 879tQ7V,8a3 

887 KRnet&t889) ETtH 
GOTO 690 

888 KRITE(6,8S0> KST,ICW 

890 RETURN , 

880 F0RMAT(//5>X»* CONGRATULATIONS t YOU DIO IT 

889 FOftMAT(//i>X,* STEPSlZfc TOO SHALL , YOU ARE 

891 FPRHATI//(SXf6D19.1I) ) 

692 F0RHAT(5Xt9£19.7/6X,3119//) 

END 


S217J 
IN TROUbLE 


• , 2020 . 12 ) 
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SUBROUTINE : DVDPF2 


CALL STATEIVIENT : DVDPF2 (Same parameters as DVDPFl) 

DISCUSSION : This subroutine is a slightly modified version of DVDPFl which 

allows usage of the output subroutine OUTGAJ. The main differences 

% 

are: 

A. Additional matrices are defined to be used as calling parameters 
for OUTGAJ. 

B. The DVOUT common area is extended to be used in OUTGAJ. 

C. The ALL common area is defined for use in OUTGAJ. 

D. An integer N is defined and read in to identify the number of 
non-zero optical rays in the next bundle of rays. 

SUBROUTINES REQUIRED : 

O, S.U. Project Library: 


1. 

DVDQ 

(JPL routine) 

2. 

FUN 

(representing FUNPL6) 

3. 

OUT 

(representing OUTGAJ) 


REFERENCES : 

Papo, Haim B. (1971). "Optimal Selenodetic Control, " 

The Ohio State University, Department of Geodetic Science , 
Report No. 156. 
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o n n 


DVDPF2 


SUBROUTINE DVDPF2tNq ,Y.YPtKU,EPtRATltRAT2 ,HHIN»HMAXtO£LT 
P tKST, lHL»PAR»FUN,OUTfKQTY:ifOTfJaH) 

SUBROUTINE FOR INTEGRATING OIFERENTIAL EQUATIONS BY THE OVDQ 
NUMERICAL INTEGRATOR . THE NUMBER OF INTEGRANDS SHOULD BE AT 
LEAST SIX 

implicit REAL*8 (A-H.O-Z) 

REAL»4 EPINQ ) ,HH IN ,HHAX t t«X tRATl ,RAT2 

DIMENSION Y( NQJ,YP( NQItYNI NC)»OT(20, NQItKDl NO)»KQ( NO), PARIS) 
P,EH3{44,4A),Al(‘«4,7S),AlH3t7S,4A.) , A3 ( 4A, 3 ) , A2 144,2) ,N{ 31) ,SEV(6) 
CCMMCfI/DVOUT/ETCB,KM»ICW,N ,ISKIP/ALL/SEV,RCM,NUT,1A,IN2,IM, I15K 
REAO(5,61)N 
N2=IN2 
KK=1H 

C KH IS THE NUMBER OF BUNDLES TO SE PROCESSED 

1SK1P=II5K 
HXSTP=SOOOOO 
KO(l)=l 
ET=PARU) 

ETFN=PAR12) 

H=PARI3) 

1CW=0 

WRITEI6,891)Y,ET,ETFN,H,0£LT 

WRI TE { 6, 892) HHIN, hMAX jRATl ,RAT2 ,NQ ,KD ( U ,MXSTP 

KLU=-l 

GOTO 881 

879 CALL OVDQ t NO ,ET , Y, YP »KO, EP , 1 FLtH,HHIN,HHAX, DELT,ETFN,HXSTP , 
PKSTjK MX,EMX,KG,YN,DT) 

GOTO 684 

883 CALL DVDQl tNQ,ET,Y,YP ,KD ,EP , 1 FL, H, HM IN ,HHAX ,OELT,ETFN,HXSTP, 
PKSTjK HX,EHX,K«,YN,DT) 

884 GOTO (88l,882»fi8ST&Si>,885,8e6,8e7,888) , IFL 

885 CALL FUN (ET,Y,YP) 

CALL OUT (ET,Y,N2,6H3,AUA1M3,A3,AZ) 

N2=NnCH + l)^2 
IF(JOa.NE.l) GOTO 683 
OELT=ETOB-ET 

IFdCW.EQ.KH) DELT=DELT + .1D0 
ICW=1CH+1 

IF(NUT.EQ.l) GOTO 888 
GOTO 863 

886 EP(KHX)=32.*EHXW£P(KHX) 

IF tCPIKHXI.GT.O. ) EP(KHX)=-£PIKHX) 

GOTO 083 
881 00 871 1=1,6 

EP(I) = Y{I)=»RATl 

IFIABSIEPm ).LT.1.E-15) EP { 1 )=-l.E-10 

871 IMEPID.GT.O.) EP(I)=-EPm 
IFINO-6 ) 874,874,872 

872 DO 873 1=7 ,NQ 
EP(I)=Y{I)*RAT2 

IFIABStEPd) ).LT. 1.6-15) 6PU)=-1.E-10 
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DVDPF2 (Cont) 


873 1F(EP(1) ♦GT.0.1 EP(1)=-EP{I) 

874 CONTINUE 
KLU=KLU+l 

882 CAUL FUN (£T,Y,YP> 

IF(KLU> 879,879t6E3 

887 HRITE{6t689) ET,H 
GOTO 890 

888 WRITEl6t880J KST,1CW 

890 RETURN 

61 F0RHAT(25I3) 

880 F0RHAT(//5Xt • CONGRATULATIONS , YOU DID IT S2I7) 

889 FGRMAT(//5X,' STcPSIIE TOO SMALL r YOU ARE IN TROUBLE ‘,2020.12) 

891 FCRHAT(//(5X,6Ul9.1in 

892 FDRMAT(5X,4E19.7/3X,3I19//) 

END 
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SUBROUTINE : EKHARD 


CALL STATEMENT : EKHARD (ET, YTH, PHILA) 

. SUBROUTINE PURPOSE : Evaluation of the Eulerian angles of the moon and the 
physical librations of the moon Including their time rates for a given 
epoch. 

INPUT PARAMETERS ; 

A. Scaler (extended precision). 

1. ET is the epoch in Julian days. 

OUTPUT PARAMETERS ; 

A. Vectors (extended precision). 

1. YTH is a 6-element vector containing the Eulerian 
angles of the moon and their time rates in the following 
order (Papo’s notation): (o, if), 6, <o, j|), 0 in radians 
and radians per day. 

2. PHILA is a 3-element vector containing the physical 
libration angles of the moon in the following order 
(Papo's notation): r, a, p in radians. 

PROGRAM DESCRIPTION : Physical libration theory as developed by Eckhardt 

fEckhardt, 19701 is used. The basic purpose of the subroutine is achieved 
by the following steps. 

A. The subroutine MEANAN.is used to compute angle .3 of the mean ■ 
orbits of the sun and the moon for the epoch ET. Then the sine/ 
cosine combination of the Delauney variables are evaluated. 

B. Eckhardt's coefficients as given in f Eckhardt, 19701 are used to 
evaluate the physical libration angles and time rates of the angles 

. at epoch ET. 

C. Combinations of the appropriate angles computed above are used 
to obtain the Eulerian angles of the moon. 

SUBROUTINES REQUIRED : 

O. S. U. Project Libi'ary: 

1. MEANAN 


REFERENCES : 

A. Eckhardt, D. H. (1970). "Lunar Libration Tables, " The Moon , 
Vol. 1, No. 2, February. 

B. Papo, Haim B. (1971). "Optimal Selenodetic Control, " 

The Ohio State University, Department of Geodetic Science , 
Report No. 156. 


Preceding page blank 
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noon 


SUBROUTINE EKHARD (ET,YTH,PHLAl 
C CALCULATION OF EULER ANGLES USING ECKHAROT'S CONSTANTS 
1HPLIC1TREAL*8 (A-H,a-Z) 

ET - EPOCH IN JULIAN DAYS 

YTH - EULERIAN ANGLES FOR THE HOON AND THEIR TI«€ DERIVATIVES 
PHI , PSI , TETA 

PHLA - PHYSICAL LIBRATION ANGLES IN LONGITUDE j NODE AND INCLINATION 
OIHENSION PHLAf3)t ' ANM(&) ,ONM C6) ,TAU{ 9 ) i CTAU(9 ) ,DTAU( 9 

PJ tSIG{5).CSlGt5)fDSIG(5) ,ROOC5) rCRQOIS) ,DRD0(5) tYTH(6) 

DATA CTAU,CSIG,CROO/1.7*9i.6,-l.A,A.2,-3.5,-16.9,1.0»15.3,lC.O,-3. 
P0,-lO.6»-23. 8, 2. 5 ,-100.6, -3. 1,-10.5,23.8, -1.9 ,-98.4/ 

P, RA0SE,P1/ 206264, 806247C96355D0, 3. 14159 

PZ653589793D0/,EIH0/, 02676900/ 

CALL HEANAN {ET,ANH,0.';H J 
C OELAUr^EY ANGLES AND THEIR TIKE RATES 
ANKH=ANH ( 2 ) -ANM { 4 1 
0NHH=DNH{2)-DNH{4) 

ANMS«=ANM{1)-ANH{3) 

DMHS=DNK Cl )-CNH {3 I 
ARLA=ANH(2)-AN'1(5) 

DRLA=ONH(2J-ONMt5) 

EL0N=ANMt2)-ANM(l» 

0L0N=DNH(2)-DNM(1) 

SLM=DSIN(ANMH) 

CLH=OCOSUNKM) 

SLS=DSIN(ANHS) 

CLS=0C05(ANH5) 

SF=DSIN(ARLA) 

CF=0C0S(ARLA) 

SD=OSINCELON) 

CD=0C0S(ELCN1 

C SINES AND COSINES OF DELAUNEY ANGLE COMBINATIONS 


S1=2.00»SF*CF 2F 

C1=CF*CF -SF*SF 

Ol=2.D0*0RLA 

S2=SLH*C1“S1*CLM t-“2F 

C2=ClH*Cl+SLH’i‘Sl 

D2=Dr«H-Dl 

S3=SLH’*C2+S2^CLM 2L-2F 

C3=CLM*C2-SLM*S2 

D3S0NHM+D2 

S4=SLH*CD-S0*CLM C-D 

C4=CLH»CD+SLH*SD 

D4=DNHH-0LQN 

S5=2.00*S4*C4 2L-20 

C5=C4*C4-S4*S4 

05=2.D0*D4 


S6-S5*CLH-SLH=»C5 

C6=C5*CLM+SLH*S5 

D6=05-DNKM 

57=S4»CLS-SLS*C4 
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EKHARD fCont> 


C 


c 


C7=CA»CLS+SLS*S4 

07=DA-DN«S 

S8=S5*C3-S3«C5 

C6=C3*C5+S3*S5 

08=D5-D3 

S9=S5*CLS-SLS*C5 

C9sC5<'CLS-«-SLS*S5 

09=D5-0NMS 

CREATION OF VECTORS FOR PHYSICAL 

TAUC1)=S8 

TAU(2)=SLS 

TAU(3)=S7 

TAUtAJ=S6 

TAU(5)=S4 

TAU (6J=SLH 

TAU(7)=S9 

TAU(8)=S3 

TAU(9)=S5 

DTAUtl)=C8 t08 

01AU(2)=CLS*DNHS 

OTAU(3)=C7 *07 

0TAU(A)=C6 *06 

0TAU(5)=CA *04 

0TAUl6)=CLf*DNMM. 


. 2F-2D 


2L-L*~20 


L18RATI0N ANGLE AND RATE CALCULATIONS 


DTAU(7)=C9 *09 

0TAU«8)=C3 *03 

DTAU(9)=C5 *05 

5IGU)=S8 

SIG(2)=SX 

SIG(3)=S2 

SIG(4)=S6 

SlGt5)=SLn 

DSIG(1)=0TAU(1) 

0SIG{2)=Cl»01 

DSIG(3>=C2*D2 

DS16<4)=DTAU<4) 

DS1G(5)=0TAU(6) 

Root 1)=C8 

RG0(2)=C1 

ROOt31=C2 

R00(4)=C6 

R0Qt5)=CLM 

DR00(1)=-S8*38 

DROO<2}=-Sl«-Ol 

DR00(3)=-S2*02 

DR00(4)=-S6*06 

DR00t5)=-SL«*nN*.H 

CALCULATION OF PHYSICAL LIBRATIUN ANGLES { TAU SICKA ROD I 

ATAU=0,D0 

ASIG^O.OO 

AROO=O.DO 
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EKHARD (Cont) 


0ATAU=0»D0 
OASIGsO.OO 
DAROO^O.OO 
DO 7I=lf9 

ATAU=ATAU+-TAU<I)*CTAUm 

7 DATAU=0ATAU+0TAU(1)*CTAUCI) 

DO 81-ltS 

ASIG=ASIG+SIG{1)*CSIG(I) 

DASIG*DASIG+DS1G(I1*CS1G( I) 

AROO*AROO+ROO ( I ) »CROO ( I ) 

8 OAROQ=OAROC+OROO(I)=*CROOtI) 

ASlG=(ASIG/EIMa J/RADSE 

DAS1G=(DAS1G/E1HQ )/RADSE 

ATAU*ATAU/RADSE 

0ATAU=DA7ALJ/RA0SE 

AROO=AROD/RADSE 

DARDO=OAROO/RADSE 

PHLA CU=ATAU*RADSE 
PHLA(2 )=ASI6»RADSE 

C CALCULATION OF EULER ANGLES (FI PSl TETA ) AND THEIR TIME DERIVATIVES 

YTHt2)=ANH(5)+ASIG 
YTH{5)=0NH(5)+DASIG 
YTH*1)=P1+ANM{2)+ATAU-YTH(2) 

YTH ( 4) =ONH ( 2 ) +DAT AU-Y TH ( 5 ) 

YTHI3)=E1MQ+AR00 
YTH{6)=0AR00 
DO 9I=lt3 
CALL REOPI (YTH(D) 

9 CONTINUE 
RETURN 
END 
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SUBROUTINE ; EPHITL 

CALL STATEMENT : EPHITL (ET, N, Q) 


* SUBROUTINE PURPOSE : The subroutine is designed to interpolate the 18 

tabulated quantities of the simulated ephemeris and to add the result 
of the interpolation to the reference case. 


INPUT PARAMETERS : 

A. Scalars 

1 . 

B. Integer 

1 . 


(extended precision). 

ET is the epoch in Julian days minus 2440000. 0 for 
which the ephemerides are required. 

(single precision), 

N is a computer unit identifier. The computer unit 
is the device on which the tabulated ephemerides are 
stored. For example, for //FT04F001 DD 
N = 4. 


C. Common Block Input (EPHEM), extended precision. 

1. ETE is a scaler which contains the epoch in Julian 
days minus 2440000. 0 of the first input data. 

2. FH is a 3-dimensional matrix (3 x 18 X 15) which 
contains ephemeris data (including second and fourth 
differences) for 14 half-day intervals beginning with 
ETE, Note that when initially calling EPHITL, 

EPHEM should contain an arbitrary set of ETE, FH 
data (normally the first block in the tabulated data. ) 

D. Data block input (extended precision) . 

1. ENCKE and PM are 7- and 6- element vectors respec- 
tively which contain constants necessary to perform trans- 
formation of the tabulated ephemeris data from perturba- 
tions to the reference case into the ephemeris output. 

AMM is a 3 X 3 orthogonal transformation matrix also used 
in the transformation. These data block variables are 
discussed further in the FUNEPH subroutine. 


OUTPUT PARAMETERS : 


A. Vector 

1 . 


(extended precision). 

Q is a vector of 18 parameters which are interpolab 
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to the desired input epoch. The variables 
contained in Q are as follows: 

- a. The state vector of the lunar mass center 


with respect to the geocentric inertially 
oriented system in km andkm/day, (six 
quantities). 

b. The Eulerian angles of the moon and their 
time rates in radians and radians/day, 
to, 0, 6 and <p, 0, 8, respectively (six 
quantities). 

•c. The Eulerian atigles of the earth and their 
time rates in the same units as above, 
denoted as 17, X, c and 77, X, I f Papo, 1971] 
Chapter 4. 


PROGRAM DESCRIPTION : The data contained on unit N is interpolated for the 
epoch input using second and fourth differences of the tabulated quantities 
and is added to the reference cases of the translatory and rotational motioi 
of the moon and the earth. 


SUBROUTINES REQUIRED : 

O, S.U, Project Library: 

1. EVERAL 

2. STVITR 

3. RE DPI 


REFERENCES : 

Papo, Haim B. (1971). "Optimal Selenode tic Control, " 
The Ohio State University, Department of Geodetic Science , 
Report No, 158. 
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EPHTTI. 


SUBROUTINE EPHITL (ET,N,ai. 

C SUBROUTINE FOR INTERPOLATING THE SIMULATEO EPHEMERIS 

C ET - EPOCH FOR WHICH THE EPHEHERIS IS REQUIRED 

C N - UNIT ON WHICH THE EPHEMEPaS RESIDES 

C Q ~ THE INTERPELLATED VALUES FOR (1 - 6) - GEOCENTRIC POSITION 

C AND VELOCITY OF THE MCON (7 - 121 — EULERIAN ANGLES OF THE MOON . 

C INCLUDING TIME RATES tl3 - 18 ) - EULERIAN ANGLES OF THE EARTH 

C INCLUDING TIME RATES 

IHPLICITREAL=»S (A-H,0-21 ^ , 

OIHCNSION Yaei,FHt3,ie,I5),ENCK.E(7),AHH(3t3) .PHIoltOda) ,STV(6) 
COMMON /EPHEH/ £TE,FH 

DATAFNCKE/ 0.220if,’000000000000 03 » 0-3796009972332398D 01 

P r 0.823611873A910790D-01 « 0.23000000000000000 00 t 

P 0.486068666695^9500 01 r 0.409160000000COOOO 00 » 

P 0.6300386L9aC00CC0D Cl/tAHH/ —0.14358831573581850 00 « 

P .0.98577089680309050 00 s -0.67396422125166200-01 t 

P 0.98961000620343970 CO j -0.14236496470816290 00 « 

p -q. 2010603008 2937690—01 t 0.73777507590927590—02 « 

P -0.B9375364E3Y11970D-01 . 0.99597068880265990 00 /. 

PPM/ 0.99999999999999990—12 t O.22O5OO0OOOC000QO0 03 t 

P 0.38979829938949710 06 , 0.54899999999999990-01' » 

P 0.24513347715837980 01 * 0.22997150219515500 00 / 

GOTO 27 

20 REAU(N) ETE»FH 
27 DIF=ET-ETE 

IFIOIF) 21,22,23 

21 rewind N 
GOTO 20 

23 IF(OIF.LT.7.00) GOTO 22 
lST=I0INT(DlF/7.00)-l 
IFdST-EQ.O) GOTO 20 

DO 24 1=1, 1ST 

24 READ IN) ETE 
GOTO 20 

22 LP=1D1NT(DIF*2. 001+1 
S5=2.DO*OMOOtOIF, .5Dul 

CALL EVERAL (SS ,LP,FH, Y , 16 , 3, 15 1 
CALL STVITR ( ET,STV,AMH,PH) 

DO 1 1=1,6 
Q{ii=Ym+sTvm 

Qll+61=Y(I+61 
1 Q(I+12)=Y(I+12) 

TT=ET-ENCKE(11 

CI71=Q(7)+EHCKE 12 ) +ENCKE (4} *TT 
CALL REDPKQt?) ) 

Q(8)=CI(8)+EMCKI (31 
Qtl0)=G{10}+fcNCKE(4) 

C I13l=Q( 13)+L.',Ci<E <5 )+ENCKE(7)»n 
CALL REDPl (Q< 1 1- 1 ) 

Q(15}=Q(15)+EriCKE(6l 

0U6)=Q116J+ENCKE(7) 

RETURN 

END 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 
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SUBROUTINE : EVERAL 


CALL STATEMENT : EVERAL (SS, LP, FDX, STV, N, Kl, K2) 

SUBROUTINE PURPOSE : Interpolation by the use of Everett's fifth order 
modified method. 


INPUT PARAMETERS : 

A. Matrices 

1 . 


B. Scalar 

1 . 

C. Integer 

1 . 


(extended precision). 

EDX(K1, N, K2), where Kl, N, K2 subscripts are 
integers (single precision); 

a. The first integer, Kl, represents the 
number of quantities needed for the inter- 
polation of each function, i, e. , the function 
value and its second and fourth modified 
differences (see Program Description). Kl 
must be set to 3. 

b. The second integer, N, indicates the number 
of functions to be interpolated. 

c. The third integer, K2, indicates the number 
of layers (see Program Description contained 
in FDX). 

(extended precision). 

SS is the interpolation factor (i. e. , it must be a decimal 
fraction between 0. 0 and 1. 0). 

(single precision). 

LP is an integer denoting the sequential number of the 
layer (out of the total of K2), preceding the point at 
which interpolation is to be performed. 


OUTPUT PARAMETERS : 

A. Vector (extended precision): 

1. STV is an N dimension vector containing the interpolated 
functions. 


PROGRAM PESO RIPTION ; Following the notation of Brouwer and Clemence on 
pages 144 and 145 TBrouwer and Clemence, 19611, the interpolation 
of value fn is given bj' (see text for notations): 

fn = fo + n6^ + Eo6^ + Ei6? E^6^ -h Ef6i + * • ♦ 


Preeeding page blank 
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To create this function the following identities are formed: 


Computer Notation Variable 

FS(1) n 

FS(2) E® 

FS(3) Et 

FP(1) n-1 

FP(2) -E| 

FP(4) -Ej 


The level is then found from LS = LP+1 and 


Computer Notation Variable 

FDX(1, N, LS) 

FDX(2, N, LS) 6? 

FDX(3,N,LS) 6t 

FDX(1,N, LP) fp 

FDX(2,N,LP) 6o 

FDX(3. N. LP) 6^ 


The interpolated variable is then {in computer format using variable 
notation) : 


f n = fin (n-1) 

+ E!6?-(-E§)6f 
+ Et6t-(-Eo)6^ 

which reduces to the Brouwer and Clemence formula. 

SUBROUTINES BEQUIRED : None 
REFERENCES : 

A. Brouwer, D. and Clemence, G. (1961). '*Methods of Celestial 
Mechanics," Academic Press, New York. 

B. O'Handley, Douglas A. et. al. , (1969). "JPL Development 
Ephemeris Number 69, " Technical Report 33-1465. 
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EVERAL 


SUBROUTINE EVERAL (SS ,LP,‘'F0X»STV,N,KItK2) 

EVERET FIFTH ORDER HOOIFIED INTERPOLATION (JPU 

FOX - COORDINATES TO faE INTERPOLATED INCLUDING SECOND AND FORTH 
MODIFIED DIFFERENCES 

FIRST SUBSCRIPT 1 - F 2 - 0EL2 3 - DEL4 
SECOND SUBSCRIPT QUANTITIES TU tsE INTERPOLATED ( N OF THEM J 
THIRD subscript SEQUENTIAL DATA SET { 3 X N LAYER ) 

STV - INTERPOLATED QUANTITIES 
SS - INTERPOLATION FACTOR 
IHPL1C1TREAL*8 {A-h,0-ZJ 

DIMENSION FS(3)»FPJ3),STV(N)tFDXIK1,N,K2) 

LS=LP+1 

SM2=SS~2,D0 

SHl=SS-l.DO 

SP1=SS+1.D0 

SP2=5S+2.00 

Fsm=ss 

FS{2)=SM1*SS*SP1/6,D0 

■FS{3)=SM2^SH1»SS’>SPI*SP2/120.DD 

FPm=SMl 

FPt23=SS»SMl=*SM2/6.00 

FP (3 ) = SP1*SS*SM I=*SH2*( SM2-1 .DO ) /120.DO 

DOA01K=ltN 

STV(KI=0.D0 

D0402I=1t3 

402 SIV(K)=STV(K)+FOX(IfK,LS)*FS(I)-FOXII,K,LP )*FP(1) 

401 CONTINUE 
RETURN 
END 
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,SUBROUTINE : F ONE PH 


CALL STATEMENT : FUNEPH (ET, Y, YP) 

- SUBROUTINE PURPOSE : Computation of the time derivatives of the perturbations 
of the state vector of the moon, the Eulerian'angles and their time rates 
for the moon and the earth in the simulated earth-moon environment describeti 
in r Papo, 1971] Chapter 4. 

INPUT PARAMETERS : 

A. Scalar (extended precision). 

1. ET is the epoch for which the time derivatives are 
required in Julian days minus 2440000, 0, 

B. Vector (extended precision). 

1. Y is a vector of 18 parameters which are the "pertur- 
b.ations" to the reference case of the simulated environ- 
ment as developed in fPapo, 1971], The order of the 
parameters is the same as the Q output vector explained 
in the EPHITL subroutine description. 

C. Common block input (extended precision). 

1. MAINFUN is a common area containing: 

a. ENCKE is a vector of 7 elements containing 
in order (see Papo, Pages 131-135, 154-159); 

(1) . The standard or zero epoch of inte- 

gration in Julian days minus 2440000. ■ 

(2) . The mean longitude of the moon. 

(3) . The longitude of the node of the lunar 

orbit. 

(4) . The mean motion of the moon. 

(5) . The mean longitude of the U axis of the 

average terrestrial coordinate system. 

(6) , The mean inclination of the equatorial 

plane of the earth with respect to the 
ecliptic system, 

(7) . The mean rotational velocity of the 

earth about the W axis. 

NOTE ; (2) through (7) pertain to epoch 
ENCKE (1) and are given in radians or in 
radians per day. 


Preceding gage blank 
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b. AMM is a 3 X 3 matrix which is an orthogonal 
transformation matrix from a coordinate system 
defined by the perigee and the plane of the 
reference orbit of the moon into the ecliptic 
system. 

c. PM is a 6-element vector containing Keplerian 
orbital parameters of the moon. The elements 
are the precision required in the solution of 
Kepler’s equation, the standard epoch in Julian 
days (minus 244. 0000. 0), the major semi- 
axis (in km), the eccentricity (in radians), 

the mean anomaly at the standard epoch (in 
radians), and the mean motion (in radians per 
day). 

2. ALL is a common area containing the vector Q. 

a. Q is an 18-element vector containing the total 
magnitudes (e. g. , reference case plus pertur- 
bations) of the ephemeris. The quantities are 
in the same order and units as listed in the Q 
output vector described in the EPHTTL sub- 
routine. 

OUTPUT PARAMETERS : 

A. Vector (extended precision). 

1. YP is an 18-element vector containing the time derivatives 
of the input Y vector. It is based on the equations of 
motion for the simulated environment developed in f Papo, 
1971] Pages 138-150 and 154-159. 

PROGRAM DESCRIPTION : This subroutine is essentially the generator of the 
simulated earth-moon environment as outlined by f Papo, 1971] Pages 
119-160. 

SUBROUTINES REQUIRED : 

A. O.S.U. Project Library: 

1. STVITR 

2. REDPI 

3. ROTATE 

4. TRIM 

5. ONEM 

B. Fortran Scientific Subroutine Package: 

1. DGMTRA 

2. DGMADD 


REFERENCES : 

Papo, Haim B. (lOVlL "Optimal Selenodetic Control. " The 

Ohio State University, Department of Geodetic Science, Report No. IHG 
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FUNEPfl 


SUBROUTINE FUNEPH lETtY.YPI 

GENERATOR OF SIMULATED EPHEMERIS (MOON AMD EARTH) 
Y( 1- A) MOON STATE VECTOR GEOCENTRIC ECLIPTIC 
Y( 7-12) MOON EULERIAN ANGLES AND TIME RATES 
Y(13-18) EARTH EULERIAN ANGLES AND TIME RATES 
POSITION t X , Y , Z 

EULERIAN ANGLES : LONGITUDE t NODE > INCLINATION 

IHPLICITREAL^a <A-H,Q-2) 


DIMENSION YUB) ,YPaa),Q{lSJ,STV{6),ENCKEt7) tXYIO) ,X{3),T1(3,3 2*P 

P) tT2(3f3) fT3(3,3) ,T4(3»3) ,TH(3,3) ,THT{3,3),TE(3t3),TET(3t3) ,VX(3) i I 

PVX3(3),XCR13*3),AMM(3,3) , GS ( 3r3 ) ,OS (3,3 ) fPS (3 ,3 ) ,UVH(3 ) ,PM C 6 J 2 

CQHMCN/MAIFUN/ENCRE,AHM,PH,NUMB/ALL/Q 

DATA ZERO, ONE »CPHL, ALF, BET , GAM, EMU, GEAR ,CONH,ES,HS ,GS 3*P 

P /O. DO, 1,00 ,.892668016, .41942130-3, .6290-3, .20957S8D-3, .327 1 

P802D-2 ,35.99241768010,3012159753997540.00 ,66067.85625,536904 2 


P6 .636, 3578363. 863 ,3*0. 00, 3579114. 284,3*0.00 ,3580615. 12o/ 
NUHB=NUHB+1 

CALL STVITR {£T,STV,AHM,PM) 

DO 1 1=1,6 
Q(IJ=Ym+STV(l) 

Q(I+6)=YCl+6) 

1 Q(I+12)=Y(1+12) 

TT=ET-ENCKEtl» 

Q(71=Q(7)+ENCKE (2 )+ENCK£ (4)*TT 
CALL REDPI10(7) ) 

Q(8)=U(8)+ENCKC(3) 

Q(10)=0(10)+ENCKE (4) 

Q (13 )=Q( 13)+6NCK£ (5 ) +ENCKE (7) *1 
CALL REDPI(0(13)) 

Q(15)=Q(15)+ENCKE(6) 

Q(16)=Q(16)+ENCKE(7) 

XX=ZERO 
XQ=ZERO 
00 2 1=1,3 
X(1)=Q(I) 

xx=xx+xm*x(i) 

2 XO=XO+STV(I)*STV(I) 

X7=0NE/XX**3.50O 

X5=X7*XX 

X3=X5*XX 

CALL ROTATE (3,Q(7),T1) 

CALL RQTATE(1,-Q(9),T2) 

CALL R0TATEt3,0(6),T3) 

CALL TRIM (T1,T2,T4) 

CALL TRIM (T4,T3,TH). 
call DGMTRA (TM,TMT,3,3) 

CALL TRIM (TMT,GS,T3) 

CALL TRIM (T3,TH, PS) 

CALL ONEM (TM,X,XY2) 
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0Hl=Tl(2,l)*T2t3,2)»Q(ll J-T1(1,1)*Q(12) 
0H2»Tltlfl)*T2(2t3)’*yai»+Tl(l,2)=»Q(12I 
flH3= T2(3,31*atll)+ QUO) 

TRE1= Al.F*{CPHL=*XY2(2)*XYZ(3)*X5-0M2»0H3)+TlU,2)*T2{3T3)’»CUl)*Qt 
P12)*TlUtl)*T2{3,2)»Qan*Q(10»+Tl(2,ll*0U2)*0ClO) 
TRE2=-BeT*(CPHl,*XYZ{ l)*XY2(3)^=X!»-0Kl*CM3)+TlUtl]*T2(3,3)*QUi)*a{ 
P12)+Tl(2f 1 J*T2(3,2)*Q(11)«C(10)-T1U ,l)*Q(12)*Q(10) 

TRE3= GAH*(CPHL*XYZ{l)*XYZ(2)=»X5-CMl»OH2)+T2t3,2;«Q( 11 J*Q<12J 
YP(10)=( Tl(l»2)*T2(3,3)*TRcl+TlU,l)»T2{3,3»*TRe2)/T2(3t2J+TRH3 
YPtll)=C-Tl(li2) ♦TRE1-Tl(l,l) +TRE2)/T2C3,2) 

.YPU2»= -Tl«l,l) ♦TRE1+T1U,2) ♦TREZ 

C 

CALL ROTATE (3,Q(13),T1) 

CALL ROTATE <1,-Q(15),T2) 

CALL ROTATE (3,Q(14),T3) 

CALL TRIM CTltT2,T4) 

CALL TRIM (TA,T3,TE) 

CALL OGHTRA CTE»TET,3,3) 

DO 23 1*1t3 
DO 23 J=lt3 

23 QSU ,J)=132135.7125D0*TE(3,I)*TEC3,J) 

CALL ONEH (TEtXtUVH) 

C 

HEL=ENCKE(7)»«0NE+EHU)-T2(3,3)*Q< 17) 

TH01=TlU,2)^T2{3,3)»0(17)*L'll6)+nCl,l)*T2{3,2)*Qll7)=^HEL+Tl{2,l) 
P*Q (18) *HE L+CEAR*U VW ( 2 ) *UVH ( 3 > *X 5 

TW02=TUl,l)*T2{3,3)=!‘Q(17)*3UU)+Tl{2Tl)*T2l3t2>*QU7)*HEL-TlUil) 
P*0 US) *HE L-CE AR-UVH ( 1 ) *UVW { 3 ) -X 3 
YPU7) = (TlU»2)*TW0i+)l(ltU»lW02)/T2(2,31 
YPU8)=-T1U,1 )vTW01+Tl (1 ,2)*rw02 
YPU6)=T2(3,2)*QU7)*QUa)-T2(3,3)*YPU7) 

C 

QL=2ERO- 
DO 3 1=1,3 

3 QL=QL+YU )*(STVU J+YU)«.5DO) 

EFO=ONE-QNE/(GNE+2.bO*OL/XO)*>*'l.5DO 

XO=ONE/XQ**l.5DO 

DO 4 1=1,3 

4 VX3U)=XO*(Ym-Xm*EFO) 

C 

DO 5 1=1,3 
DO 5 J=l,3 

XCR{I,J)=~2.50C*XU)*X( J)*X7 
IF(I.NE.J) GOTO 5 
XCRU,J)=XCRU,J)+XX#X7 

5 CONTINUE 

CALL OGHAOO (PS,0S,T1,3',3 ) 

CALL TRIM (XCR,T1,T2) 

DO 6 1=1,3 


- 38 - 



FUNEPH (Cont) 


DO 6 J=lt3 

IFCI-NE.J) GOTO 6 

T2{ I ,J)=T2(I»J )+X5’«‘(ES+HSJ 

6 CONTINUE 

CALL ONEM (T2,XtVX) 

DO 7 1=1,3 

YP(I+3)=-C0NM=i'(VXm+VX3{I) J 
Y?(I )=Y(I+3 > 

YP(I+6 )=Y(I+9 ) 

7 YP{I+12)=Y(I+15) 

99 RETURN 

END 
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SUBROUTINE : FUNPL5 


CALL STATEMENT ; FUNPL5 (ET, Y, YP) 

SUBROUTINE PURPOSE : The subroutine calculates the second time deriva- 
tives of the physical libration angles of the real moon and the state 
transition matrix and the parameter sensitivity matrix used in the numer- 
ical integration of the physical librations of the moon as described in 
rPapo, 1971] Chapter 3. 32. The selenodetic positions of the earth and 
the sun are obtained from the JPL DE-69 tape TO' Handley, 1969]. 


INPUT PARAMETERS : 


A. 


B. 


Scalars (extended precision). Also note the Common Area 
Parameters section. 

1. ET is the epoch in Julian days for which the derivative? 
are needed. 

Vector (extended precision). 

1. *Y is a 60-element vector containing, in order (supplied 
by the calling program DVDPF): 

a. The physical libration angles and time rates 
(r,u,p,f,a,p). 


b. The 36-elements of the state transition maU- 
stored columnwize. 

c. The 18 elements of the parameter sensitivit 
matrix stored columnwise. 


OUTPUT PARAMETERS : 

A. Vector (extended precision). Also note the Common Area 
Parameters section. 

1. YP is a 60-element vector containing the time derivative 
of the input Y vector: 

YP 5 Y 


COMMON AREA PARAMETERS : The program extensively uses common area 
input/output, listed separately in this section. The common areas are 
MAIF, EXE, ALL, CETBLl, CETBL2, CETBL4. 

A. Area /MAIF/ALF, BET, GAM, TEQ. 

1. ALF, BET, GAM are the rations between the 
moments of inertia of the moon: 


a = 


C-B 




C-A 

B 


, V = 


B- A 
C 


Precedins page blank 
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rPapo, 19711. They are extended precision scalars 
evaluated in the main program and passed through the 
subroutine calling FUNPL5 (usually DVDPF). 

2. TEQ is the mean inclination of the lunar equator with 
respect to the mean ecliptic of date (I). 

B. /EXE/TO, B, SUN, NUMB. 

1. TO is an extended precision scalar containing the epoch 
of the previous step in the integration at which FUNPL5 
was called. 

2. B is an extended precision 3-element vector containing the 
geocentric position vector* of the moon in kilometers, 

3. SUN is a vector similar to B, containing the geocentric 
position vector of the sun at epoch TO. 

4. NUMB is a dummy integer. 

C. /ALL/SEV, A, NUT, KLU. 

1. SEV is an extended precision 6-element vector containing 
the mean longitude and longitude of node of the moon's 
orbit and their first and second time derivatives 
respectively, evaluated at the epoch ET. 

2. A is an extended precision dummy scalar. 

3. NUT is a single precision integer which is normally 

set to zero. The integer is assigned the value one and the 
integration is terminated if abnormal program operation is 
detected. 

4. KLU is a single precision integer which is used to control 
the program actions. If set to one, the partial deriva- 
tives of the YP elements 7-60 are integrated with the 
physical libration angles. If set to an integer greater 
than one, only the physical libration angles (elements 
1-6 of the Y vector) are integrated. 

D. Common areas /CETBL1/CETBL2/CETBL4/ are required for 
input from the JPL DE-69 tape [O'Handley et al, , 1969]. 

PROGRAM DESCRIPTION : The FUNPL5 program is designed to find the second 
time derivatives of the real moon's (e. g. as described by DE-69) motion 
[Papo, 1971] Chapter 3 . The routine requires the use of subroutine 
TRIK and the logic used in the program is given below. 

A. The geocentric positions of the moon and sun are read in from the 
DE-69 tape. 

B. The position vectors from (A) are transformed into selenocentric 
positions in the mean ecliptic system of date (ET). 

C. The Eulerian angles of the moon and their time derivatives 

are evaluated from combinations of the Y and SEV vector elements. 

D. The TRIK subroutine is called to evaluate the time derivatives 
of the physical libration angles and to create (in Papo's notation) 
the 0 and Q matrices need for the evaluation of U and Q 

r Papo, 1971] Page 9b. 
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REFERENCES : 

A. O' Handley, Douglas A, et al. (1969). 'tTPL Development Ephemerls 
Number 69J* Technical Report 32-1465 . 

B. Papo, Haim B. (1971). "Optimal Selenodetic Control, " The 
Ohio State University, Department of Geodetic Science Report Nu. 
156. 
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FUNPL5 


SUBROUTINE FUNPL5 (ET,YrYPJ 

SUBROUTINE FUNPl.5 FOR READING THE OE-69 TAPE , CALLING THE TRIK SULRGUTINE 
AND ADMINISTERING THE FORMATION OF DERIVATIVES OF THE PHYSICAL LIBRATION 
ANGLES AND THE PARTIAL DERIVATIVES MATRICES (TRANSITION AND SENSITIVITY) 
1HPLICITREAL*8(A-H,0-Z) 

DIMENSION Y{601*YP(60J,SEV(6) t Z ( 6) i 2P(3) »DC3) |B( 3) »SUM(3) t 

P0ME{3),F{4) ,TEl{6),TE2{fa)fIR£a3) ,S (6, 12) ♦TU3,3r»T2{ 3»3) tT3(3»3) 
P,TETA(6,6),FETA(6,3) 

Y(l-6) - PHYSICAL LIBRATIONS ANGLES 

Y(7-A2)- STATE TRANSITION MATRIX IN COLUMNS (TAU ...RODOT) 

YI43-60) - PARAMETERS SENSITIVITY IN COLUMNS (C22t BETA,C20 ) 
COHKDN/HAIF/ALFfBET »GAM,TEO /EXE/TO » B f SUN*NUM3/ALL/ SEViAtNUTt 
PKLU/CETBL1/AU,REA ,TPOf EMR/CE T5L2/ICH, ICE » 1 RE/C£TBL4/St F 
DATA PItTSE /3.141592o535897VD0t0.D0/ 

C ET ~ EPOCH IN { JO-2R4OOO0.) 

IF(ET.EQ.TO)' GOTO 2 
T0=ET 

UJ=ET+ 2440000 .DO 
CALL READE (UJiTSEflER) 

C READING FROM THE DE-69 TAPE 

IFtlER.NE.O) GOTO 6 
DO 1 1=1»3 
D(I)=S(I,11) 

1 OME(I)=S{IflO) 

C POSITIONS OF EARTH(D) AND SUN (ONE) IN MEAN EQUATORIAL OF 1950 

CALL HEANAN (UJ.TE1.TE2) 

CALL PRECSS (UJ»T2,T3,1) 

CALL ROTATE ( 1 ,7E 1 1 6 ) ,T1 ) 

CALL TRIM (T1,T2,T3) 

CALL ONEH (T3,0»r>) 

CALL ONEH (T3,0H6»SUN) 

C POSITIONS OF EARTHO) AND SUN(SUN) IN MEAN OF DATE ECLIPTIC 

SEVtl)=TElt21 
SEV{2)=TEU5) 

SEV(3)=TE2(2) 

SEV(4)=TE2(5) 

SEV 1 5) =— 0. 2965A088Q— 13+0.A059^23430—20*{ ET+24980.D0) 

SEV(M= 0.543658266D-13+0.477579229D-20*IET+24960.00) 

C SEV VECTOR OF MEAN LONGITUDE , NODE AND THEIR FIRST AND SECOND DERIVATIVS 

2 CONTINUE 

ZU)=Y{1)*“Y(2)+SEV(1)-SEV(2)+PI 
CALL REDPI (Z( U ) 

. 2(2)=Y(2)+SEVC2) 

CALL REDPI (Z(2)) 

Z(3)=Y(3)+TE0 

Z(4)=Y(4)-Y(5)+SEV(3)-SEV(4) 

Z(5)=Y(5)+Sev(A) 

Z(6)=Y(6) 

C Z - EULERIAN ANGLES CALCULATED FROM Y 

CALL TRIK (Z,2P,TETA,FETA) 

DO 3 K=l»3 
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YP{K)=Y(3+K> 

3 YP{3+K)=2P(K) 

C YP{l-6) ASSIGNMENT OF DERIVATIVES FOR THE PHYSICAL LI8RATI0N ANGLES 
IFtKLO.GT.lJ GOTO 17 

C KLU IS SET GREATER THAN I FOR THE CASE WHEN NO PARTIALS ARE GENERATED 

DO 16 1=1,6 
00 13 J=l,6 
I J=J*6+I 
YP(IJ)=O.DO 
DO 13 L=l,6 

13 YPf IJ)=YP(IJ)-»-TETA(l »U^Y<J*6+U 

C CREATION OF U DOT IN VECTOR FORM YP(7-42> 

DO 15 J=l,3 
IJ=36+I+6*J 
YP(1J)=0.D0 
DO 14 L=l,6 

14 YPIIJ)=YP{IJ)+TETA(I,L)*Y(36+6*J+L) 

15 YPC1J)=YP(1JJ+FETAU ,J) 

C CREATION OF Q OOT IN VECTOR FORM YP 143-60 J 

16 CONTINUE 

17 NUHB=NUHB+1 
RETURN 

6 HR1TE(6,70)ET,NUHB 
NUT=1 
RETURN 

70 FORM AT (5X, 'SOMETHING IS WRONG IN THE DATA', F20. 8,1101 
END 
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SUBROUTINE : FUNPL6 


CALL STATEMENT ; FUNPL6 (ET, Y, YP) 

' SUBROUTINE PURPOSE : The subroutine accomplishes the same purposes 

FUNPL5 except that the simulated environment is used in lieu of the 
JPL ephemeris DE-69. 


INPUT PARAMETERS : Input parameters are exactly the same as for FUNPL5. 

OUTPUT PARAMETERS : Output parameters are exactly the same as for FUNPL5. 

COMMON AREA PARAMETERS : Common storage is the same as FUNPL5 except 
that the areas CETBL1/CETBL2/CETBL4 are not required- 


PROGRAM DESCRIPTION : The computation performed by the subroutine is 

exactly the same as FUNPL5 except that the simulated ephemeris data is 
read from a data set No. 4 through the use of EPHITL. 


SUBROUTINES REQUIRED : 

A. O. S. U. Project Library; 

1. REDPI 

2. TRIK 

3. EPHITL 


REFERENCES : 

Papo, Haim B. (1971). "Optimal Selenodetic Control, " 
The Ohio State University, Department of Geodetic Science 
Report No. 156. 


— 
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SUBROUTINE FUNPL6 tET«Y,YPj 
IHPLICITREAL#a fA-H,0-2} 

DIMENSION Y(60) ,YPt60J,SEV{&) , Z( 6 J t 2Pf 3» tD{ 31 ,B ( 3 ) tSUN (3 ) , 

P0HE(3),F(4)tTcl(6)tTE2(6) ,1 RE tl3), St 6, 12), TlC3,3),T2I3t3),T3C3,31 
P,T£TA(6,6),FETA(6,31,0U(1GJ 
C Yll-6) - PHYSICAL HBRATIONS ANGLES 

C Y(7-A2J- STATE TRANSITION MATRIX IN COLUMNS tTAU ...RODOT) 

C Y{43-60) - PARAMETERS SENSITIVITY IN COLUMNS (C22 ,UETA ,C20 J 

C0HHCN/MAIF/ALF,BET,GAH,TEQ /EXE/TO,8,SUN,NUH8/ALL/ SEV,A,NUT, 

■ PKLU 

DATA PIiTSE /3. 1415926535897900,0. DO/ 

C ET - EPOCH IN (JD-2440000.) 

IFtET.EQ.TO) GOTO 2 
TO=ET 

CALL EPHITL tET,4,QU) 

DO 1 1=1,3 
B(I)=QU(I) 

1 SUN(I)=1.D 05 

2 CONTINUE 

Zll)=Ym-Y(2)+SEV(l)-SEV(2)+Pl + (SEVt3)-SEV|4n*tET-222.5D0) 

CALL REDPI (Ztin 

Z(2)=Y{2)+SEVt2)+SEVt4J#{ET-222.5DOJ 
CALL REDPI (Z(21) 

Z(3)=Y(3)+TECJ 

Z(4)=yt4)-Y(51+SEV(3)-SEV{4) 

2l5)=yt5)+SEV(4) 

Z(6)=Y{6) 

C Z - EULERIAN ANGLES CALCULATED FROM Y 

CALL TRIK (Z,ZP,TETA,FETA) 

DO 3 K=l,3 
YPtK)=Yt3+K) 

3 YP(3-*K)=2P(K) 

C YP(l-6) ASSIGNMENT OF DERIVATIVES FOR THE PHYSICAL LIBRATION ANGLES 
IF(KLU.GT.l) GOTO 17 

C KLU IS SET GREATER THAN 1 FOR THE CASE WHEN NO PARTIALS ARE GENERATED 

DO 16 1=1,6 
DO 13 J=l,6 
IJ=J*6-*-I 
YP(1J)=0.00 
DO 13 L=l,6 

13 YP(1J>=YP{IJ)+TETA{I,U*Y(J*6+L1 

C CREATION OF U DOT IN VECTOR FORM YP(7-42) 

DO 15 J=l,3 
IJ=36+I+6*J 
YP(1JJ=0.D0 
DO 14 L=l,6 

14 YPI IJ)=YP( IJ) +TETA( I, U*Y(36+6*J+U 

15 YP(1J)=YP(IJ)+FETA(1,3) 

C CREATION OF Q DOT IN VECTOR FORM YP( 43-60 : 

16 CONTINUE 

17 CONTINUE 
RETURN 
END 
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SUBROUTINE : FUNST3 


CALL STATEMENT : FUN STS (ET, Y, YP) 


SUBR OUTD^E PURPOSE : B3iSica.lly the subroutine accomplishes two purposes: 

A. The subroutine evaluates the time derivatives of the perturba- 
tions of the 18 elements of the simulated ephemeris (similar to 
FUNEPH). 

B. The subroutine evaluates the time derivatives of the state 
vector of a satellite in the simulated earth moon environment. 


INPUT PARAMETERS : 

A. Scalar 

1 . 

B. Vector 

1 . 


(extended precision). 

ET is the epoch in Julian days at which the time derivative 
are to be evaluated. 

(extended precision). 

Y is a 24-element vector containing, sequentially: 

a. The 18 perturbation of the elements 
of the simulated ephemeris. as 

defined by comment cards in the statement 
listings. Units of the parameter sets are 
kilometers, kilometers per day, radians 
and radians per day respectively, 

b. Elements 19-24 contain the state vector 
(kilometers and kilometers per day) of the 
satellite in the selenocentric system. 


OUTPUT PARAMETERS : 

A. Vector (extended precision), 

1. YP is a 24-element vector returned from the subroutine 
containing the time derivatives of the Y vector described in 
the above paragraph. Units are per 6a.y time derivatives 
of the original vector. 


COMMON AREA PARAMETERS : Common block storage is used to transmit infor- 
mation between the subroutine and the main program or calling subroutine. 
Two common areas are used; 

A. Area /MAIFUN/ENCKE, AMM, PM, NUMB. 

1. The extended precision variables ENCKE, AMM and 
are quantities which arc described in the FUNEPH 
Program Description, 
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2. NUMB is a single precision integer which serves as a 
counter to indicate how many times the FUNST3 sub- 
routine has been called. 

B.- Area ABL/Q. 

1. The extended precision vector Q contains the 18 elements 
of the simulated ephemeris in the same order as in the 
Y vector. To obtain the Q vector the perturbed elements ' 
of the Y vector are added to the "reference” case. 


PROGRAM DESCRIPTION : The theoretical basis of.PUNST3 formulation is outlined 
in PPapo, 1971] Chapter 4, Section 4.4. 

SUBROUTINES REQUIRED : 

A. O. S.U, Project Library: 

1. STVITR 

2. REDPI 

3. ROTATE 

4. TRIM 

5. ONEM 

B. Fortran Scientific Subroutine Package: 

1. DGMPRD 

2. DGMTRA 

3. DGMADD 

REFERENCES : 

Papo, Haim B. (1971). "Optimal Selenodetic Control, " The 
Ohio State University, Department of Geodetic Science, Report 
No. 156. 
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FUNST3 


SUBROUTINE FUNST3 tET,Y,YP) 

SIMULATED ENVIRONMENT 

SATELLITE TRAJECTORY GENERATOR COWELL EQUATIONS OF MOTION FOR SATELLI £ 

NO RANGE OR RRATE DATA GENtRATED 

YC 1“ 6) MOON STATE VECTOR GEOCENTRIC ECLIPtiC 

Y( 7-12) HCON EULERIAN ANGLES AND TIME RATES 

Ytl3-i8) EARTH EULERIAN ANGLES AND TIME RATES 

Y(19-2A) SATELLITE STaTE VECTOR SELENOCENTRIC ECLIPTIC 

POSITION : X , Y , Z 

EULERIAN ANGLES : LONGITUDE t NODE ^ INCLINATION 

IMPLICIT REAL*8 (A-H,0-Z) 


DIMENSION Y(24) ,YP ( 24 J tCM 18 ) t ST V 1 6) , ENCKE ( 7 ) , 5(3 ) ,R( 3) , }< ( 3 ) , T1 ( 3, 3 5*P 

P)»T2(3»3)tT3(3»3) ,T4( 3 t3 ) ,TH( 3,3) ,THT(3 ,3 ) ,TE (3,3 1 ,TET (3 ,3) , VX( 3) , 1 

PVX3(3),VS{3),EHS(12),XH(12,3) ,SM(12,3) ,S SM (12 ) , SC0(3 ) , 0HSU2 ) ,PM 1 6 2 

P ),RRRR{3,4 ) ,GS(3,3),QS(3,3),PS(3,3),XYZ(3),UVWt3 4 

P),XCR(3,3),SCR(3,3) ,RCR ( 3 , 3 ) , VR (3 1 , AHM (3 , 3 ) ,VR3(3) 5 

COMMON /HAIFUN/ENCKE,AHH,PH,NUMB /ALL/ 0 

DATA ZERO, ONE, EMM, CPHL, ALE, BET, GAM, EMU, GEAR, CONS, ES,HS 5*P 

P /O. DO, 1.00,81. 300, .C926680 16, .4194213D-03 , .6290-03 

P, .2095788D-03,.3278C2D-C2, 3699241 76800. 00,. 2973560 16,66067.836250 1 

PO, 5369046. 63600/ 

DATA GS,CQNM / 3578363.86300,3*0.00,3579114.26400,3*0.00,3560 

P6 15. 12600,3012159753997540.00/ 5 

DATAEHS,XH/.234556D“06,-.246487D-C6,.1891820-C6 3 

P, -.1462610-06, .14437oD-06, -. 157250-06 ,. 347270-07, -.743310-07, 4 

P. 1235040-06, -.990660-07, .998140-07, -.1025630-06, 5 

P1619.177, 1532.288,1457.013,1630.737, 938. 9‘»0 , 1 110 .463 ,750.1 23, 6 

P836. 942, 713 -644, 582. 937, 2*730. 123, -433. 527, 559-617, 666. 747, 1.0, 7 

P-940.490, -1325. 401, -1302.581,-223.928, 856.851, lo09.474, 1304. 8 

P561, 435.527, 449.355,593.968,-151.986,-594.968, -1117.734, -151. 9 

P986, 868.554, 1504.745, 1530.968, 301.319,-869.554,-1505.743/ 10 


NUHB=NUHB-»1 

CALL STVITR (ET,STV,AHH,PH) 

00 1 1 = 1,6 
Qm=Y(l)+STV(I ) 

Q(I+6)=Y(1+6) 

1 Q(I*12)=Y1I+12) 

TT=ET-ENCKE(1) 

Q(7)=Q{7)*ENCKE(2)+ENCKE(4)*TT 
CALL RE0PI(0(7)) 

CI(8)=Q(8)+ENCKL (3) 
QU0)=Q(10)->£NCKE(4) 

Q C13)=Q{13)+ENCKE (5 ) +ENCKE ( 7 ) *TT 
CALL REDPnQtl3)) 
Q(15)=0(15)+ENCKE(6) 
QU6)=0C16)+ENCK.B{7) 

DO 22 1=1,12 
22 SSHU)=ZERO 
C 

XX=ZERO 
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XD=2ER0 

SS=ZERO 

RRsZCRO 

DO 2 1*1,3 

Xll)*Q(II 

SCI>=Y{2+I8j 

RII»=X(I)*Sm 

SCQ(l»=2ERO 

xx=xx+xa}*x{ij 

XO=XC+STVU>*STVU) 

SS=SS*SCI>»S(1» 

2 RR=RR+R{1)*R(I> 

X7=0NE/XX**3 » 50 0 
S7=ONE/{EHM*SS**3.5DO) 

•R7=0NE/RR **3.500 

X5*X7»XX 

S5=S7*SS 

R5=R7*RR 

X3=X5*XX 

S3=S5*SS 

CAUL ROTATE (3,Q{7),TU 
CALL R0TATEU,-'CH9),T2) 

CALL R0TATEC3tQI8>,T3) 

CALL trim m,T2,T4) 

CALL TRIH (T4,13,TM) 

CALL OGHTRA (TH,TMT,3,3) 

CALL TRIH {THT,GS,T3) 

CALL TRIH {T3,TH,PS} 

CALL ONEH (TM.XtXYZ) 

CALL 0GMPRO (XH,TH,SH,i2,3,3) 

OHI=TU2,U*T2(5,2)*QCIU-T1{1,1J*0(121 

0H2=Tia,n*T2(2,3)*0aU*Tm,2)*Q<12) 

0M3= T2C3,3)*Q{lXl-»’ Q{10) 

TREl* ALF»<CPHL*XYZl2>*XYZt3)*X5-OM2*OH3J-«-Tl{l,2)*T2(3t3)*i3t 1U*Q( 
P12l+Tm,U*T2n,2)*a(lll*Qa0}+Tl(2,lJ*Qfl2)*Qfl0) 
TR£2®-8£T*(CPHL»XYZm*XYZ<3)*X5~0Hl*0.M3}+TI(l,l)*T2{3,3»*QCllJ*a( 
PX2»*Tlf2,l}*T2l3,2}*oai>*CJtlO)-Tm,l)*C{12}*atlf') 

TRE3* GAM*(CPHL*XY21 1)*XYZ{ 2) «X5-0M1*0H2 >+T2 (3,2 )»Qni J*Q ( J2) 
yPUOl=C TlU,2l*T2{3,3>*TREl+TZa,U*T2f3,3)*TRE2}/T2(3,23+TRE3 
YPUI) = (-T1 IIt 2J *TRE1-T1(1,U *IBE2J/T2t3,2) 

YP(12)= -TlU.l) *TRE1+TI(1,2) *TRE2 

CALL ROTATE {3,Q(13J,T1) 

CALL ROTATE a,-0(15),T2 
CALL R(3TATE (3,ca4),T3> 

CALL TRIH {T1,T2,T41 
CALL TRIM {T4,T3,TE) 

CALL OGHTRA (TE,TET,3,3J 
DO 23 1*1,3 
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FUNST3 (Cont) 


C 


C 


C 


00 23 J=l,3 

23 OS { I , J 1=1321 35 . 71 25D0»TE (3,I1*TE13,J) 
CALL ONEM (TEjXtUVH) 


HEL=ENCKE(7)»(CNE+EMU)-T2 (3 ,31*0(17) 

THQ1=T1( 1,21 *T2{ 3,3 1*0(17 1*0(18 l+Tl ( 1, 1)*T2(3 
P*Q 1 1 8 1 *H£L+CE AR*U VW ( 2 1 *U VW ( 3 1 *X 5 
TK02=T1(1»1 1*72(3,31*0(17 1*0 (18)+T1( 2,1) *T2 (3 

P*0( 181 *HEL-CEAR*UVW ( 1 1*UVW t 31 »X5 
YPtl7)=(Tl(l ,2)*THC/1+Tl(l»i 1*TW021/T2(2,31 
YP(lfcl=-Tl(l,ll*1W01+Tl(l,21*lW02 
YP( 161=T2(3,21*0( 171*Q( 18 1-12 (3,3 1*YP (17 1 


,2)*qa71*HEL+Tl(2, 

,21*0(17)*HEL-1XU» 


II 

11 


QL=7ER0 
OR=ZERO ■' 

00 3 1=1,3 
DO 10 J=l,12 
SH(J,I)=S(I1-SH(J,I1 
10 SSH(J1=SSH(J1*SM(J, 11**2 

QR=QR+S(I 1*(X( I ) + S( 1 1*0.5001 

3 QL=QL+Y{I1*(STV(1)+Y(I)*.500) 
EPR=(3NE-DNE/( ONE+2.DO*QR/XX)*’ 
EFO=ONE-ONE/(r;tI6 + 2-DO*QL/XO)*' 
XO=ONE/XO** 1.500 

DO 4 1=1,3 

VR3{I)=X3*(S(1)-€FR*R(I1 1 

4 VX3(I)=X0*(Y(I )-X(I)*EFO) 

DO 15 1=1,12 

OHS(I)=EHS( I )/S5H(l 1**1.500 
DO 15 J=l,3 

15 SCO(Jl=SCU(Jl+DrfS(I)*SH(I,Jl 


DO 5 1=1,3 
DO 5 J=l,3 

XCRU, J)=-2.5D0*X(11*X( J)*X7 
SCR(I,J1=-2.5D0-?S(I)*S{J)*S7 
RCRd, Jl=-2.5UO*R( 11*K( J)*R7 

IFd.NE.Jl GOru 5 
XCRd,J)=XCRd ,J1+XX*X7 
SCRd,J)=SCR(I.ai+SS*S7 
RCR d , J 1 =RCR ( I , J 1 +RR*R7 
5 CONTINUE • 

CALL DGHADD (PS, OS, 71,3, 3) 
CALL TRIM (XCR, 71,72) 

CALL TRIH (SCK,PS,T3) 

CALL TRIH (RCR, OS, 74) 

DO 4 1=1,3 

DO 6 J=l,3 

IFd.NE.J) GOTO 6 

72d,J 1=72(1, J)+X5*(ES+HS1 

73(1 , J 1=73 (I , J 1 +S5*HS*S3 
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FUNST3 (Cont) 


T4(I,J)=T4(I»J)+R5*£S 


6 


7 

99 


CONTINUE 

CALL ONEM {T2,X,VX» 

CALL ONEM (T3*S,VS) 

CALL ONEM (TAtRtVR) 

00 7 1=1*3 

YPI I+3)=“C0NM«^(VX(1)+VX3( I 
yP(I+21)=-C0NS*{VR(I)+VStI 

YPd )=Y(I+3 ) 

YPd+6 )=Y(I+9 ) 
YP(I+12)=Y(I+15) 
YP{I+ia)=Y(I+21) 

RETURN 

Ef4D 


)) 

J+SC0(I)-VX(I)+VR3(I 


)} 
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SUBROUTINE : GASTIM 


CALL STATEMENT : GASTIM (JD, DOBLQ, DLONG, GAST) 

SUBROUTINE PURPOSE : The subroutine computes the Greenwich apparent 
sidereal time at a given epoch, given the nufation in obliquity and 
longitude for the epoch. 

INPUT PARAMETERS : 

A. Scalars (extended precision). 

1. JD is the epoch in Julian days minus 2437000. 0. 

2. DOBLQ is the nutation in obliquity in radians. 

3. DLONG is the nutation in longitude in radians. 

OUTPUT PARAMETERS : 

A. Scalars (extended precision). 

1. GAST is the Greenwich apparent siderial time for the 
epoch specified by JD. 

PROGRAM DESCRIPTION : [Fajemirokun, 1971] Page 51. 

SUBROUTINES REQUIRED : 

O. S.U. Project Library: 

1. MEANAN 

REFERENCES : 

Fajemirokun, F. A, (1971). "Application of New Observational 
Systems for Selenodetic Control, " The Ohio State University, 
Department of Geodetic Science, Report No, 157 . 
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uu u 


GASTIM 


SUBROUTINE GASTIHI JOf OOBLQiOLONGtGAST) 
lHPLICITP.EAL*fatA-H,0-I) 

DOUBLE PRECISION JO 

DIMENSION ANM(6),0NM(6) ,,,nn/ 

DATA PI2/6. 233 18S307l79366n0/fRAD0EG/57. 2957795130823200/ 

THIS ROUTINE COMPUTES THE G.A.S.T. 

INPUT PARAM. — JOtNUTATIomiN RADIANS) IN LONG £ 0§LQ 
OUTPUT IS CAST IN RADIANS 
FJOs'JO'*-2437COO.ODO 
C ALLHE ANAN { F JO t ANH, DNH) 

08LQ=ANH(6) 

OBL0=O6LQ+D0BLQ 


EQE=OLONG<'DCOS { QB LQ ) 
IJO=IDINT( JD) 
DJO=DFLOAT(IJD) <0.500 


FRAC=JO-DJD 

1F{FRAC.LT,0.D0)FRAC=1.D0+FRAC 


UT=FRAC*2A.OO 

TH=IFJD-2A150Z0. 001/36525. 00 

GMST=UT+ (6*D0+36 • 00/ 60^ D0+A5-B36D0/5600*D0) ♦( C6A0164* 
♦1 *TH+( 0 .092900/3600.00) *TH*TH 
GHST=0M0D t GH5T , 2A .DO 1 
GHS7=GHST*15.00/RADD£G 
GAST=GHST+EQE 
GAST=DH0DtGAST,PI2) 

RETURN 


54200/3600.00 


60 F0RHATI//»10X»025.161 


ENO 
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SUBROUTINE : GEULAN 


CALL STATEMENT ; GEULAN (JD, THETA, PHI, PSI, DOBLQ, DLONG) 

SUBROUTINE PURPOSE : Computes the three Eulerian angles (0, >J), <o) between 
the mean ecliptic system of 1950. 0 and the average terrestrial system. 

INPUT PARAMETERS : 

A. Scalars (extended precision). 

1. JD is the Julian date minus 2437000. 0. 

2. DOBLQ is the nutation in obliquity in radians. 

3. DLONG is the nutation in longitude in radians. 

OUTPUT PARAMETERS ; 

A. Scalars (extended precision). 

1. THETA is the Eulerian angle 0 in radians. 

2. PHI is the Eulerian angle (O in radians. 

3. PSI is the Eulerian angle j|) in radians. 

PROGRAM DESCRIPTION : The development of the earth's Eulerian angles is 
outlined in Report 157, Section 2.35. 


SUBROUTINES REQUIRED 


O. s.u. 

Project Library; 

1. 

PRECES 

2. 

NUTATE 

3. 

GASTIM 

4. 

STRANS 

5. 

MEANAN 

6. 

ROTATE 

7. 

TRIM 


REFERENCES : 

Fajemirokun, F. A. (1971). "Application of New Observational 
Systems for Selenodetic Control, " The Ohio State University, 
Department of Geodetic Science, Report No. 157. 



noono 


GEULAN 


SUBROUTINE GEULANl JD,TH£TA»PHIt PSl tOOBLQ,DLONG} 

IMPL1CITR£AL« (A-H,0-Z) 

OOU6LEPRECISIONJD 

DIMENSION PRE{6) ,PRECH(3,3),0NUTA{3»3J,SHAT{3,3JtANM(6UDNM(6|, 
«REll3,3)tTH3,3),T2{3,3),A(3,3J 
data ET/2A33282,A23D0/ 

DATA PI2/6.28318550717958600/ 

THIS ROUTINE COMPUTES THE 3 EULERIAN ANGLES BETWEEN THE MEAN ECLIP 
TIC SYSTEM U950.0) AND THE MEAN TERRESTIAL ROTATING SYSTEM 
EXTERNAL ROUTINES REOD. -PRECESt NUTATE iGASTIHtSTftANStMEANAN, ROTATE* 
TRIM. INPUT PARAH. — JUL. OATE,NUTAT.IN OSLO. AND LONGITUDE 

F JD= JD+2437000 .ODO 
CALLPRECES(FJO»PRE,PRECMJ 
C ALLMJTATE ( F JD * D08L0 1 OLONG, ONUTA ) 

CALLGASTIH(JOtOQBLQfOLONG,GAST) 

CALLSTRANSJFJO.GASTjSHAT) 

CALLMEANANl ET,A^JM ,ONM ) 

OBLQ-ANM(6) 

CALLROTATEd.-OSLCtREl) 

CALLTRIM(PRECM*R£1*T1) 

CALLTR1H(DNUTA,T1 ,T2) 

CALLTRIKCSMAT,T2,A) 

THETA=DARCOStA(3,3) } 

$T=OSIN( THETA J 
PSI=DARC0S(A(3,2)/ST) 

P$IsPI2-PSI 

PHI=UARCOSJ (Al l,l)*A(3»2)-All,2)*AC3tlJ )/ST) 
DS=-(IA(2tl)*A(3*2)-A(2*2)»A(3,l) )/ST) 

IF (DS*LT.O, DO) PHI=P 12-PHI 
50 FORHATC/, (3025.16)) 

60 FORMAT!//, 2X.5D20. 12) 

RETURN 

END 
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SUBROUTINE : IN VS PE 


CALL-STATEMENT : INVSPE (A, B, M, N, K, L) 

, SUBROUTINE PURPOSE : Copies a layer of a three dimensional input matrix (A) 
into a two dimensional matrix (B). 

INPUT PARAMETERS : 

A. Scalar integers (single precision). 

1. M, N, K are dimensions of the input matrix A 
(rows, columns, layers). 

2. L is the layer of A (consisting of K layers) which is to be 
copied into matrix B for output. 

B. Matrices (extended precision). 

1, A is the M X N X K input matrix. 

OUTPUT PARAMETERS : 

A. Matrices (extended precision). 

1. B is the output MXN matrix containing layer L 
of the input matrix A. . 

PROGRAM DESCRIPTION : None 

SUBROUTINES REQUIRED : None 

REFERENCES : None 
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mvsPE 


c 

SUBROUTINEINVSPE{ A,8»M,N»K,L) 
C MAKING B EQUAL TO LAYER OF A 

IMPLICITREAL=«‘a (A-H,0-Z) 
D1HENSI0NA{H*N,K) 

DOli=lfM 

002J-lfN 

Z B(I,JI=Aa,J,LJ 
1 CONTINUE 
RETURN 
END 
C 
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SUBROUTINE : LADIS 


CALL STATEMENT ; LADIS (ET, PE, PM, XE, XM, XCE, X, D, ZD) 

SUBROUTINE PURPOSE : The program computes the distance between an 
earth observatory and a lunar reflector at a given epoch and the 
zenith distance of the ray at the earth observatory. The program 
presupposes the observations are made in the simulated environment 
described in fPapo, 1971] in a refraction free environment assumed in 
[Fajemirokun, 1971], 


INPUT PARAMETERS : 

A. Scalar (extended precision). 

1. ET is the epoch of observation in Julian days. 

B. Vectors (extended precision, 3 element). 

1. XE is the position vector of the earth observatory 
in a body fixed system (km) . 

2. XM is the position vector of the lunar reflector 
in a body fixed system (km) , 

3. XCE is the geocentric position of the moon in the 
simulated environment (km), 

C. Matrices (extended precision). 

1, PE is a 3X3 matrix used to transform coordinates 
from the earth fixed system to the mean ecliptic 
coordinate system. 

2, PM is a 3 X3 matrix used to transform coordinates 
from the moon fixed system to the mean ecliptic system. 


OUTPUT PARAMETERS : 


A. 

Vector 


1. 

B. 

Scalars 


1. 


2. 


(extended precision). 

X is the topocentric position vector of the lunar 
reflector (km). 

(extended precision), 

D is the distance between the earth observing station and 
the lunar reflector (km). The variable is set to zero 
if the zenith distance precludes observation. 

ZD is the zenith distance of the lunar point in degrees. 


PROGRAM DESCRIPTION; The program logic follows the formulation given in 
r Fajemirokun, 1971] for laser ranging. 
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A. First, the topocentric position vector of the lunar reflector 
is obtained. 

B. The simulated laser range is found from the topocentric vector. 

C. The zenith distance of the observed range is computed. 

D. If the observing altitude of the ray is between 30 and 70 degrees, 
control is returned to the calling program. If the "observation’' 
cannot be made, the distance is set to zero and control is 
returned to the calling subroutine. 


SUBROUTINES REQUIRED : 

A. O.S.U. Project Library: 

1. ONEM 

B. O.S.U. Utility Library: 

1. MADD 

2. MSUBT 


REFERENCES : 

A. Fajemirokun, F. A. (1971). "Application of New Observational 
Systems for Selenodetic Control, " The Ohio State University, 
Department of Geodetic Science, Report No. 157. 

B. Papo, Haim B. (1971). "Optimal Selenodetic Control, " The 
Ohio State University, Department of Geodetic Science. Report 
No. 156. 
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LAPIS 


C 

SUBROUTINE LADIS{tT,P£,PM,Xe,XM,XCE,X,D,ZO) 

C THIS PRL.GRAH SIMULATES RANGES 10 1H£ KClON 

C PROGRAM USES SIMULATED ENVIRG.\Ht.\T DATA (CREATED 6Y H. PAPOJ 

C £1 IS EPOCH OF CbSERVAHON IN JULIAN DAY 

C P£ MATRIX transforms FROM DV.T TO MEAN ECLP. SYSTEM 

C PH matrix TRANSFORMS XYZIKCON) TO MEAN ECLP. SYSTEM 

C ■ XE,XM — POS- VECT. OF EARTH PT. AND MOON PT. IN UVW C XYZ SYSTEM 

C XCE — GEOCENTRIC (ECLP.) POSITION OF SELENOCENTER 

C O IS THE COMPUTED DISTANCE 

C X — ^THE rOPOCENTRlC POS. OF MOON POINT IN MEAN ECLP. SYSTEM 

C ZO IS THE zenith DISTANCE CF MOON POINT 
IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION PE(3,3),PM(3,3) tX£(3),XH(3),XCE(3 ) tXEE13) ,XH£ (3) ,XEH(3) » 
»X(3) 

DATA RADDEG/57. 29577951306232/ 

DS1=DSIN(20.D0/?ADU£G) 

OS2=OS1N(60-DO/RAODEG) 

DS3=0CDS (80.DO/RADDEG) 

c 

CALLCNEH(PE,X£,XEE) 

CALLDNEH(PHtXH,XME) 

C ALLKADD ( XCLtXHE , 3 , I »XEM) 

CALLHSU6T(XEHtXE£t3.1»X) 

D2=X(I)*X(1)+X(2)*X(2)+X13)*X{3) 

D=0SQRT(D2) 

C the following CHECKS IF OBSERVING ALTITUDE IS 30<ALT.<70 

A3=0.00 
EB=0.00 
A=0.00 
B=O.DO 
E=0.00 
D050l=1.3 
AD=Aft+XEE{ 1 )*X( 1 ) 

EB=EB+XME ( U* (-X ( I ) ) 

A=A-fXEE(i)*XEE( 1) 

B=B+Xin*Xt I) 

E=E+XME(n*XME( I ) 

5C CONTINUE 
A=DSQRT(A) 

6=DS0RT{B) 

E=DSQRT(E) 

CST=AB/(A*b) 

CH=6B/(E»fa) 

«D=OARCUS(CST) 

2D=0D=»RADDEO 
IF(CST) 60,60,65 
65 CONTINUE 

IF(CST,L£.DS1.0R.CST.GE.DS2)GOT060 
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T.ADIS fCont) 


THE FOLLOWING. CHECKS IF MOON PCiNT IS ..-O DEG. OFF LUNAP. LlHb 
1 F(CH.LT.DS3)GDT060 
GCT0999 
oO CONUNUe 

3 impossible to ObSERN 

D=O.DO 
99v RETURN 
Vo F0RMAT<3D25.12) 

END 
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SUBROUTINE : LASOLV 

CALL STATEMENT : LASOLV (IDP, XP, IDL, XL, B, EMINV, QLA, 

PX, TEMP, LVEC, MVEC, EN, COREL, NE, NM, NU, NOB, NP, 

NL, BM, W) 

SUBROUTINE PURPOSE : This subroutine is the basic routine for adjusting 
simulated laser ranging. 

INPUT PARAMETERS : 

A. Integers {single precision) from main program. 

1. ' NE is the number of earth station coordinates. 

2. NM is the number of lunar reflector coordinates. 

3. NOBS is the number of laser observations. 

4. NP is the number of earth stations. 

5. NL is the number of lunar reflectors. 

B. Integers (single precision), card input within subroutine. 

1. IDP is an NP element vector containing the 
earth station numbers. 

2. IDL is an NL element vector containing the lunar 
Station numbers, 

3, IDl is an input earth station number which an 
observation is made from. 

4, ID2 is an input lunar station number to which an 
observation is made, 

C. Matrices (extended precision) card input -within subroutine. 

1. XP is an NPX3 element array containing the geocentric 
positions of the observing stations in kilometers. 

2. XL is an NLx 3 element array containing the seleno- 
centric positions of the laser retroflectors in kilometers. 

D. Scalars (extended precision) card input within subroutine. 

1. ET is the epoch in Julian days minus 24370000.0 days 
of an observation. 

2. X is the topocentric position vector of the lunar 
reflector at ET (in kilometers). 

3. D is the simulated lunar range in kilometers at ET. 

4. ZD is the zenith distance in degrees of the simulated 
observation at ET. 

E. Scalars (extended precision) input from disk within subprogram. 

1. ETE is the initial epoch for lunar data (see EPHITL), 
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2. ETl is an input epoch for earth Eulerian angle data. 

3. ET2 is an input epoch for lunar Eulerian angle data. 

F. Vectors (extended precision) input from disk within subprogram. 

1. ET and El are element vectors described in. the - 
PLADIS Program Description. 

2. Y is a 60-element vector described in the FUNPL5 
Program Description. 

G. Matrices (extended precision) input from disk within subprogram. 

1. FH is a 3 X 18 X 15 array described in the EPHITL 
Program Description. 

2. SE is a 3 X 6 array described in the PLADIS 
Program Description. 

output parameters : Many of the extended precision variables in the 

Call Statement' are included for Fortran object time dimensioning. The 
important output from the routine is threefold: 

A. The array QLA dimensioned Nux NU is printed and represents 
the variance-covariance matrix of the adjustment. The variances 
are listed in TFajemirokun, 19713 Page 200. 

B. The array COREL dimensioned NU xNU is the correlation 
matrix resulting from the adjustment F Fajemirokun, 19713 
Pages 202-204. 

C. The W vector dimensioned NU is the solutior vector obtained 
from the adjustment 

PROGRAM DESCRIPTION * Tbe program carries out the computations outlined 
theoretically in [Fajemirokun, 19713 Chapters. 


SUBROUTINES REQUIRED : 

A. O. S.U, Project Library: 

1. EPHITL 

2. PMAT 

3. LADIS 

4. PLADIS 

5. MATPA 

B. Fortran Scientific Subroutine Package: 

1. DMINV 

2. DGMPRD 

C. O. S. U. Utility Library: 

1. MWRITE 

2. MSCALE 

3. MDUP 

REFERENCES : 

Fajemirokun, F. A. (1971). "Application of New Observational 
Systems for Selenodetic Control, " The Ohio State University, 
Department of Geodetic Science. Report No. 157. 
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LASOLV 


subrouting LASCLVI IDPtXP» lDLtXL*8»EHINV ,QtA»PX»TEHPtLVEC»KVEC*£N» 
♦COREl.,Nf'tNNiNU,NOB»NPiM.»6H»W) 

IMPLICIT REAL*e«A-H,U-Z) 

DIMENSION I DP (NP J » XP (NP 1 3 )jIDL(NL)»XLlNLt3) » SJNOBtNU) »EHINV (NQo) t 
>»QLA{NU»NU ) t PX(NU ) »TfcMP ( NOS ) » LVEC ( NU) »HVcC( NU) t LNt NUtNU) * 
*COREUNU,NUI»FH(3,lo,lS),CU(lo> ,XCF< 3) ,xm3) .XE13) »X (3 J ,PE( 3r3) , 
*PK(3,3)fEPt6),EL(6},S£t3,6)»SH(3.6),PS«(3,3) t£>l(21)tS2(3J,Y{60), 

*Z1(6) ,Z2(6) .E(15) 

DIMENSION BM{NUfNCSj,W(NOBJ 
COMMON/EPHEH/ETctFH 
1 CONTINUE 

CON=206.264b062 470964 DO 

ETC=O.DO 

00251=1. NE 

PX(I)=i-00/t25.DO*25.DOJ 

25 CONTINUE 
Kl=NE+7 
KZ=N£+6+NM 
D021I=KI,K2 

PX( I )=1. 00/(1000. 00*1000 -DO) 

21 CONTINUE 
K3=K2+1 
K4=K2+3 
00261 =K3,K4 
Jl=I-NM-6 

PX(Jl)=l-DO/( 1.00*1.00) 

PXIJl+3) =1.00/ (0.500*0. 500) 

PXtn=1.00/( 20.00*20.00) 

PX( I+3)=l.00/( 10.00*10.00) 

26 CONTINUE 
K5=K4+3 

PX(K5+l)=l. 00/(0. 500*0.500) 

PX(K5+2)=1.DO/(2.0DO*2.CDO) 

PXIK5+3)=1. 00/ (0.0100*0.0106) 

DU29I=1.3 
0028 J=l, 3 
PSH( l»J)=0.u0 
21! COfUINUc 
DC29K=1.6 
f.E(I»K)=O.DO 
SM(I.K)=0.00 
29 CONTINUE 
D0301=1,NP 

READ(5.el)I0P(I),(XPU.J) .J=l»3) 

0027 J= 1.3 

1VALUE=1DINT(XP( I,J)*10.D0) 

XPUiJ)=DPLOaT( IVALUE) 

XP(1.J)=XP(I,JI/10.00 
27 continue 

WRITE (o, 81 )iOP( 1) .(XPll.JJ »J=1.3) 
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LASOLV (Cont) 


30 continue 
0040I=ltNL 

READ(S,&i>l!jL(Ut(XLntJ) tJ=lr3) 
0035J=1,3 

lVAt.UE=IDlNT(XH I, JI ) 

35 XLU»J)*OFLOAT( IVALUE) 

WRITE(6»&IH0LU>7(XL(ltJ}>J=lt3) 

• 40 CONTINUE 
REWIND 4 
REWIND 2 
REWIND 3 
REA0t4) ETE,FH 
IJ«10 
11=16 
NR=0 
L1=NE*6 
L2=NE+6+NM 
L3=L2+6 

45 READ(5»92)ET,lDltID2,X,01tZ0 
IFCET.EO-O.DOJOOTU70 
lF(ET.6Q.ETC)G0TU5e> 
CALLEPHITUET,4,UU) 

00501 =li3 
XCE(l)=OUm 
50 CONTINUE 

57 READ(31 ETltFPf { IStU»K) tX=l»61 tl=li31 
IFtCTl.NE.ETlGCTOS-T 

CONTINUE 

THETAE=EPm 

PSI6=EPI2) 

PH1£=EP(3) 

CALI.PHAIIPS1E»THETA£,PHIE,PE> 

58 REA0(2)ET2tEE»Y 
JFIETZ.NE.ETIGOrOSS 
CONTINUE 
00531=li3 
0053J=1»6 
JJ=J»6+I 

53 SHa#J)=Y(JJI 
00541 =i» 3 
D054J=1,3 
JJ=J*6-«-I+36 

54 PSHII.J»=YUJ) 

THETAH=EU1) 

PS1M=EL(2) 

PHlM=eH3) 

CAU.PMATIPSlH,lHETAH,PhIH,PM) 

56 CONTINUE 
D060J=lt3 
XHU1=XL(ID2,J) 

XetJ)=XPUOl,J) 

60 CONTINUE 
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CALLLA01S( tT,PttPM,XEtXM,XCE,X,D»ZO) 

WRITE (6f9A)ET, 101, 1D2 ,X,D,ZD 

IF(O.EU.O,DO)GOT0999 

NR=NR+1 

W(NR)={01-D)*10Q0.00 

CALLPLADlS(ET,Xfc,XH,XCfc»X,Dl,Pt,PH»EP,EL,Se,SM,PSH,Bl,a2J 

J=ClDl*3)-2 

K=(I02*3)-2 

B(NR,J)=-B1(1) 

B(NR,J+1)=-61(2) 

b(NR,J<-2l=-bH3J 

B(NR,K+Ll)=-Bi(10) 

6(NR,K+Ll+i)=-BHll) 

B(NR,K+Ll+2)=-Bl(12J 

00651=1,6 

B(NR,I-fNE)=-eiCI+3)/CQN 
B(NR,1+L2)=-&U I+12J/C0N 
65 CONTIAUt 

BINR,L3+1»=-61( 19)/C0N 
BtNR,L3+2)=-BU20)/UGN 
B(NR,L3+3)=-81(21J/C0« 

ETC=6T 
GQT0A5 
70 CLNT1\UE 

HR1TE(6,95}NR 
CALLHHRITElW, l,.'JOfa-, »0W •» 

CALLHATPA(6,NR,NL;,tH J.WtOLA, TtHP) 

DO 75 1=1, NU 

QIA (1 , 1 )=CLAU , I J +PX U ) 

75 CCKTINUE 
FACTCR=0.00 
00641=1, NU 

64 FaCTOR=FACTCR*DLCG10 { CLA( 1 , I ) ) 

FACTOR=1.DO/10.0G**( f ACI0R/DFL0AT<NU> ) 

CALLKSCALE IFACrUR,t,l.A,NU,NU) 

CALLKDUP(OLA,NU,NJ,£N) 

CALLO«I!<V{OLA,NU,OfcI ,LVeC,HVeC) 

CALLDGHPROC EN,(.LA,COr<.E:LtNU,.NJ,NU) 

CAULHSCALE C FACTOR ,CiLA ,.NUtMJ ) 
CALLHWRITE(C0R£L,NU,NU,*1CHKM 
00621=1, NU 
D062J=l,NU 

62 C0REL(1,J)=QLA( 1 , J ) / » DSCP. T( Qt.A( 1 , 1 ) J *DSORT {gLA ( J , J) ) ) 
CALUH WR I TE 1 COREL , NU , NO , * 1 CC.1 • ) 

WRIT£(6,90) 

N=0 

0076i=l,NE,3 

N=N+1 

VAl=QLA(lyl) 

VA2=<.LAU+l,I + l) 

VA3=QLAH+2,I+2) 
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76 WRlTE|6,96)NtVAl,VA2tVA3 
WRIT£C6t90) 

M=0 

Ll=U+X 

0Q77I=Lltl-2.3 

N=N+1 

VAi=CLA(lf 1) 

VA2=(JLAU + ltI+l) 

VA3*0LACl+2,l+2) 

77 WRITE (0,96) N, VAX, VA2.VA3 
KR1TE(6,90) 

Ll=Ll-l 

00781=1,6 

feU)=QLAiI+NE,I+NU 

ECI+6)=UI.AU+L2,H-l2) 

78 CONTINUE 
E(13)=QLAIL3>l,L3+l) 

E { 19 ) =OLA ( L3+2, L3+2 ) 

E U5 ) =OLA ( L3+3, L3+3 ) 

WR1T£(6,97)E 

WR1TE(6,90) 

WRITE (6, 120) ( (CCiRELl I, J),J=1,9) ,1=1,9 ) 
HR1TE(6,120)UCCR£L( I , J 1 , J=10, 10) , 1 = 10, 18) 
HR1TE(6,120M (CCREL11,J),J=19,27),I=19,27) 
WR1TE(6,90) 

HRITE(6,120)nC0F.EL(I,J),J=l,9),I=10,18) 
WRlTE{6,l20)((ajRfcUl,JJ,J=l,9>,l=19,27) 
HRl TE ( 6, 120 H (COREU I, JJ,J=10, 18 1,1=19,27) 
MRITE(6,90) 

WRITE(6,200)( (uLAU, J) ,J=1,9),I=1,9) 
MRITE(6,200)UQLA(1,J),J=10,18),1=10,18) 
WRITE (6,200 )((QLA( I, J),J=19, 27), 1=19, 27) 
WR1TE(6,90) 

WR1TE(6,200) ( (QLAU, J) ,J=i,9 ) , 1=10,18) 
HRITE(6,200)( lOLA 1 1, J),J=1,V 1,1=19,27) 

WRI TE(6, 200 ) C (Ui.A( I, J), 1=10,18), 1=19,27) 
HR1T£(6,90) 

0C73J=l,NU 

D0731=1,NR 

73 BH(J,1)=8( 1,J)^EHINV( n 

CAU.0GMPR0(6«,W,EmINV,NU,NR,1 ) 
CALLDGHPRDtQLA,EHlNV,H,NU,NU,l) 

CALLHSCALE (-1 .00,rt,NU, 1) 

HR1T£(6,2C2)W 

WR1TE(6,90) 

HR1TE16,99)DET 

WR1TE(6,96) 

999 RETURN 
B1 FORMAT (15,3020.9) 

62 F0RHAT(I2,aX,3F10.9) 

90 FORMAT (IHl) 
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92 F0RMAT(F7.3,2I2,4F1o.7,F5«1) 

94 FORHAT(/,S>X»F7«2»2I5,‘tD19.iO,5X,F5«l) 

95 FORMATC//,I10,//) 

96 F0RMAT(15X» I5,5X,3Di7.7) 

97 F0RMAI(20X,3D17.7) 

9b FORiMAK//, 5 X 7 'CONGRATULATIONS, JOb IS OVER •) 
99 FORMA T(//,D25.i4) 

100 F0RHAT(//(5X,12F7.3) ) 

HO F0RMAT(//(/,5X,6f-7.5) } 

120 FORMAT! //( 5X,9F7.3) ) 

200 F0RMAT(//<5X,ll>9Di0.2} > 

202 FORMAT <lHl,///(/,SX,6Dl4. 6) ) 

END 
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SUBROUTINE : LUC A 


CALL STATEMENT : LUCA (N, D, DAN) 


SUBROUTINE PURPOSE : Creates special skew matrices for differentiation 

of a rotation matrix as outlined by Lucas (t is an independent variable). 


dKN(0!) 

dt 


Rn(a) I'M 


da. 

dt 


T T, /V da 

Rn (a) 


INPUT PARAMETERS : 

A. Integer N denotes the rotation axis (N = 1, 2, 3). 

r?rv 

B. Scalar D is the derivative of the rotation angle . 

' ‘ T? //V \ 

If D = 1. 0, the DAN matrix will yield — ^ 

• • •• = . ha 


OUTPUT PARAMETERS : 


A. 


DAN is 
D (i. e. 


a 3X3 matrix of partial derivatives multiplied b 

ha 


PROGRAM DESCRIPTION : 

A. The DAN matrix is zeroed. 

B. Indicies Nl, N2 are set as (by integer authentic) : 

N Nl N2 

1 .2 3 

2 3 1 

3 1 2 

C. The following elements of the DAN matrix are set as: 

DAN (Nl, N2) ' r D 
DAN ( N2, Nl) = -D 

Resulting, for axis 1, in: 


Preceding page blank 




DAN 


0 

0 

0 


0 

0 

-D 


0 

D 

0 


For axis 2 : 


- DAN 


0 0 -D 
0 0 0 
D 0 0 


For axis 3 : 


DAN 



D 

0 

0 


0 

0 

0 


SUBROUTINES REQUIRED ; None 


REFERENCES ; 

Lucas, James (1963), "Differentiation of Orientation Matrix, 
Photogprammetric Engineering, July. 


LtJCA 


SUBROUTINE LUCA <N,0,DAN) 
REAL*8D,DAN(3,3) 

0011=1,3 
001J=1 ,3 
1 0AN{ I, J)=0,00 

N1=N+1-(<N-H)/A)#3 

N2=N+2-({N+2)/A-)*3 

0AN(M1,N2)=0 

DAN(N2,Ni)=-D 

RETURN 

END 
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SUBROUTINE ; MATPA 


CALL STATEMENT : MATPA 


(A, lA, JA, P, ANS, TEMP) 


SPB ROUTINE purpose : The subroutine forms the product A'PA rvhere A is 
an lAx JA matrix and P is an lA vector representing a diagonal lAxIA 

matrix and the product is dimensioned JAXJA. 


INPUT PARAMETERS : 


A. 

Matrix 


. 1. 

B. 

Vector 


1. 

C. 

Integers 


1 . 

2 . 


(extended precision). 

A is the input matrix which is dimensioned lA x JA. 
(extended precision). 

P is a vector with lA elements which represent the 
diagonal elements of an lAXIA matrix. 

(single precision). 

lA is the number of columns in the A matrix. 

JJ is the number of rows in the input A matrix. 


OUTPUT PARAMETERS ; 

A. Matrix (extended precision). 

1. ANS is an JAx JA matrix representing the product A' PA. 

B, Vector (extended precision). 

1. TEMP is a work vector contained in the I/O parameters 
for object time dimensioning, 

P ROGRAM DESCRIPTION ; Denoting the ANS matrix by N the subroutine computes: 


N = A'PA 


SUBROUTINES REQUIRED ; None 
REFERENCES ; None 
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SUBROUTINE MATPAC A, lA, JAvP,ANS,TEMP) 

DOUBLE PRECISION A( i A, JA) »P <1 A) .ANSI J At JA) , TEMPI iA) 

DO ^0 J=1.JA 

DU 25 I=i,lA 

TEMPI 1)=0.0D 00 

TEMPI 1)=PI 1)*AM.J) 

25 CONTINUE 

DO 40 1=1. J A 
ANSI J.I)=0.00 00 
DO 40 K=1»1A 

ANS(J.I)=ANS(J»I)+AIK.I I’i^TEMPIK) 

^0 CONTINUE 
RETURN 
END 
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SUBROUTINE : MCROSS 

CALL STATEMENT ; MCROSS (X, Y, R) 

SUBROUTINE PURPOSE : The routine computes the.cross product R of two input 
vectors X and Y. 


INPUT PARAMETERS ; 

A. Vectors (extended precision). 

1, X is a 3-element input vector. 

2. Y is a 3-element input vector. 


OUTPUT PARAMETERS : 

A. Vector (extended precision). 

1. R is the cross product of the input vectors (3 elements). 


PROGRAM DESCRIPTION : The routine computes R = XXY 


SUBROUTINES REQUIRED : None 


REFERENCES : None 
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MCRQSS 


SU3RCUTINE HCROSS {X , Y » R ) v r v 

COMPUTES THE CROSS PRODUCT — R — OF TWO MATRIX VECTORS X S Y 

IMPLICIT REAL=«'8(A-H,0--Z) 

DIM£NSIONX(3),Yt3) ,R(3) 

RU)=tXC2)’>Y(3))“(X(3)=i‘Y(2) ) 

RC2)=tX(3)=>=Y(l) )“(X(1)’^Y(3) > 

R{3)=CXC1)=^Y{2) )-(X(2)=«‘Y( 1) ) 

RETURN 

END 
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SUBROUTINE : MEAN AN 


CALL STATEMENT ; ME AN AN {ET,. ANM, DNM) 

SUBROUTINE PURPOSE : Calculates elements of the mean orbits of the moon and 
the sun and their time derivatives for a given epoch. 

INPUT PARAMETERS : 

A. ET is the epoch in Julian days. 

OUTPUT PARAMETERS : 

A. ANM is a 6-element vector (extended precision). All elements 
are in radians. 

1. ANM(l) is the longitude of the sun. 

2. ANM(2) is the longitude of the moon measured from the 

mean equinox along the ecliptic up to the ascending node 
of the lunar orbit and then along the orbit. 

3. ANM(3) is the longitude of perigee of the sun. 

4. ANM(4) is the longitude of perigee of the moon measured 

along the ecliptic and the lunar orbit as ANM(2). 

5. ANM(5) is the longitude of the ascending node of the 
mean lunar orbit. 

6. ANM(6) is the obliquity of the mean equator of the earth 
with respect to the ecliptic, 

B. DNM is a 6 -element vector (extended precision). All elements 
are in radians per day. Each element in DNM is the time deriva- 
tive of the respective element in ANM. 

PROGRAM DESCRIPTION : 

A. The number of Julian days since Jan. 0, 5, 1900 is computed 
from: 

TD = ET - 2415020.0 

B. The number of Julian centuries since Jan. 0. 5, 1900 is computed 
from: 


TC = TD/36525. 0 = TD/DINC 
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C. The polynomial expressions for T^, T^, are formed in 
conjunction with 1/57. 295* • • = 1/RADD as follows: 

T(l) = l./RADD 
T(2) = T/RADD 

T(3) = t’/RADD 

T(4) = T^/RADD 

D, The SSP matrix multiple routines are used to form the mean 
angle in ANM and their time derivatives in DNM: 


sANMi 

sDNMi 


eClVIi 4TX 
sDC 4 - 4X1* 


_1 

36525 


E. The ANM vector is examined using the subroutine REDPI to 
reduce all angles to the interval (0-2ir)., 

F, The DNM vector is multiplied by the. time derivative T as below: 

DNM - DNM ( 30252 X 24 X 60 X 0. 07436574^ 


SUBROUTINES REQUIRED : 

A. Fortran Scientific Subroutine Package: 

1. DGMPRD 

B. O. S. U. Project Library: 

1. REDPI 

REFERENCES : 

Mendez, J. C. and R. J, Stern (1969). "Geographic and 
Selenodetic Coordinate Transformation Program" TRW Note 
No. 69-FMT-749, Project Apollo Task, NSC/TRW A-193. 
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n o o o n o o 


MEANAN 


SUBROUTINE MEANAN (ET,ANM,DNM) 

CALCULATES MEAN ANGLES {RAOIANS) AND TKEIR TIKE DERIVATIVES 
(RADIANS PER EPHEHERIS DAY) ; THE ANGLES ARE MEASURED FROM THE 
MEAN EQUINOX OF DATE (ST) 

ET - EPOCH IN JULIAN DAYS 

THE CONSTANTS UF CM MATRIX GENERATE MEAN ANGLES AS FOLLOWS: 

SUN LCNGf MOON LONG, SUN PERIGEE LONG, MOON PERIGEE LONG, 

MOON NODE LONG, EARTH EQUATOR OBLIQUITY 
IHPLICITREAL*8(A-H,0-Z) 

DIMENSION CM(6,A) ,CC (6,A ) ,T(A),ANH(b) ,CNM<fe) 

DATA CM /279. 696677800, 270. ^3A1 63900, 28 1-220R33300,33A.3295556DO 

P, 259, I 83275D0, 23. A5229AAOO, 36000. 76692500, A81267.6S3IA1700,1-. 7 1917 
P5000, 4069. 034033300,-1934. 142008300, “.C1301 2500, .30250-03, -.113333 
P3D-C2, .45277780-03, -.01032500, ,20777730-02 16388890-05, .000,. 188 
pa 90-05,. 3333330-05, -. 125D-0A,. 222220-05,. 5027778D-C6/,0C/360O0. 768 


P92500, 48 1267.883 1417D0, 1.71917500, 4069,03 A0333D0,-1V3A^. 142006300,- 6 

P.013012500, .6050-03 ,-.2266670-02,. 9055560-03, -.0206500,. 41555560-0 7 

P2, -.32777780-05 , .000, .566670-05,1,0-05,-, 3750-04, .666670-05*. 15083 8 ■ 

P30-05,6*0.DO/,OINC /36525.D0 /,RA0D/57. 295779 9 

P5 13082320900/ 10 


TD=rET -2415020.00 
TC=TD/DINC 
T(1)=1.D0/RA00 
T(2)=TC*T(1) 

T(3)=TC*T(2) 

T(4)=TC*T(3) 

CALLDGHPRD(CH,T,ANH,6,4,I ) 
CALL0GHPR0(0C,T,DNM,6,4,1) 
0011 = 1,6 

CALL REOP I (ANMm) 

1 DNHtI)=DNM(I)/OINC 
RETURN 
END 
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SUBROUTINE : NUTATE 


CALL STATEMENT ; NUTATE (JD, DOBLQ, DLONG, DNUTA) 

SUBROUTINE PURPOSE : Computes the nutation matrix for transformation from 
the mean to the true celestial Cartesian coordinate system of date. 


INPUT PARAMETERS 


A. 


Scalars 

1 . 

2 . 

3, 


(extended precision). 

JD is the epoch in Julian days. 

DOBLQ is the nutation in obliquity in radians. 
DLONG is the nutation in longitude in radians. 


OUTPUT PARAMETERS ; 

A. Matrices (extended precision). 

1. DNUTA is a 3 X3 vector containing the nutation matrix. 


PROGRAM DESCRIPTION : The nutation matrix N is formed from 

N = Ri(-€-A€) R3(-M) Ri(€) 

where c is the mean obliquity (obtained from the subroutine MEANAN), 
Ac is the nutation in obliquity and A0 is the nutation in longitude. 

SUBROUTINES REQUIRED : 

O.S.U. Project Library: 

1. MEANAN 

2. TRIM 


REFERENCES : 

Mueller, Ivan I. (1969). "Spherical and Practical Astronomy 
as Applied to Geodesy, " Frederick Ungar Publishing Co. , 
Now York. 
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non 


NUTATE 


SUBROUTINE NUTATE (JD»UObLQ,OLGNGtONUTA) 
lNPLICITREAL*fa tA-H,D-Z) 

DOUBLE PRECISIDN JD 

DIrtEMSION ANM{6) tU'NM(6) ,RTE113,3) iREl (3,3 ) ,RDPSI (3,3),Tl(3t3), 
=<'DNUTA(3,3) 

THIS ROUTINE CDHP. THE NUTATION MATRIX FOR TRANSF. FRH MEAN TO IRuE 
REQO. INPUT PARAM. — JULIAN DATE, NUTATION IN LONG AND QBLw. 
ROUTINES USED ARE HEANAN,TR1H 
CALLHEANAN( JD,Ai'JM,DNH) 

0BLQ=ANH{6) 

TOBLQ?OaLQ+DOBLQ 
C ALLROTATE Cl,-TOB L0,RTE1) 

CALLROTATE(l,OBLQ,REl) 

CALLROTATE(3,~DLONG,RDPSI ) 

CALLTRIM{RDPSI,REl,Ti ) 

CALLTRIM(RTE1,T1,DNUTA) 

RETURN 

60 FCRMAT{/,5X,4D25.16) 

END 
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SUBROUTINE : ONEM 


CALL STATEMENT ; ONEM (A, B. C) 


SUBROUTINE PURPOSE : Multiplies a matrix A by a vector B resulting in a 
vector C. 

INPUT PARAMETERS : 

A. MATRIX A (3X3) 

B. VECTOR B (3) 

OUTPUT PARAMETERS : 

A. VECTOR C (3) 

PROGRAM DESCRIPTION : 

A. aC^ = 3 A 3 3 B 1 

SUBROUTINES REQUIRED : None 

REFERENCES : None 
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ONEM 


t. 

SUBRCUTINE ONEM <A,B,CJ 

C PRODUCT OF MATRIX iA) AND VECTOR (6) RESULTS IN VECTOR (C) 

R£AL=5'8A<3,3) ,B(3) ,C(3) 

CCl)=AUtlJ’^B{l)+All,2)-B (2)+A{ I,3>=5'3(3) 
C(2)-A{2,i>’^B<l)+A(2,2>=»=B(2)+A12,3)-BC3J 
C C3) =A(3,1)=S'B{1I+A(3,2)*9 (2)+A(3,3)=i'8 (3) ' 

RETURN 
■ END 
C 
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SUBROUTINE : OPTOBS 


CALL STATEMENT ; OPTOBS (ET, X, JOK, APP, BUNDLE) . 

SUBROUTINE PURPOSE : The subroutine generates a bundle of rays simulating 
optical observations from a point in space exterior to the lunar surface 
(e.g. , the earth or a satellite) to an array of SO points defined in 
fpapo, 1971] Pages 162-168. 

INPUT PARAMETERS : 

A, Scalars. 

1. ET is an extended precision variable which contains 
the Julian date at which the observations are to be 
simulated. 

2 . JOK is a single precision integer indicating the number 
of the earth observatory from which the optical observa- 
tion was made (JOK = 1, 2, 3) or JOK = 9 for satellite 
based observations. 

3. APP is an extended precision variable which contains 
the sine of one-half of the field angle of the camera used. 

B, Vector (extended precision). 

1. X is a 9-element vector which represents: 

a. The geocentric position of the moon or the 
selenocentric position. of the satellite (as 
applicable). 

b. The Eulerian angles in radians of the seno- 
centric system (elements 4 - 6 ). 

c. The Eulerian angles in radians of the mean 
terrestrial system (elements 7-9). All 

X quantities are referred to the inertially 
oriented coordinate systems. 

OUTPUT PARAMETERS : 

A. Matrix (extended precision). 

1. BUNDLE is a 2 X 30 matrix which contains 30 angular 
observations in radians of lunar points from a projection 
center. The first index indicates the angle u and the 
second index indicates the angle ju referred to the 
Bi, Bg, B 3 reference system as described in fPapo, 1971] 
Pages 164-165 and Figure 4. 7. 

COMMON AREA PARAMETERS : 

A. Area /OPTO/TB, WTER, SSUN, ETER, PCSP, STI, ST2, 

STS, IFLAG. 
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1. TB is an extended precision 3X3 array supplied 

by the main or calling program and is used to transform 
coordinates from the simulated inertial system to the 
Bi, Ba, Bs system defined in rpapo, 1971], 

2. WTER is an extended precision variable evaluated in 
the subprogram and is the selenodetic longitude of the 
terminator west of the subsolar point in degrees. 

3. SSUN is an extended precisoxn variable evaluated in the 
subprogram and represents the selenodetic longitude of the 
subsolar point in degrees. 

4. ETER is an extended precision variable evaluated in 
the subprogram and is the selenodetic longitude of the 
terminator east of the subsolar point in degrees. 

5. STl is an extended precision variable which is supplied 
by the main program and is, used to check if a particular 
point is at least 5 degrees from the terminator and on the 
lighted side of the disk of the moon. 

6. ST2 is a variable similar to STl and is used to check if the 
point on the moon as seen from the projection center is 

at least 20 from the limb of the moon. 

7. ST3 is a variable similar to STl and ST2 and is used to 
check if the lunar point is at least 20° above the station 
horizon (in the case of earth based observations). 

8. IFLAG is a single precision integer used as a flag by 
the subroutine. The flag is set to one if the observing 
station is in daylight and no observations are made. 

The flag is set to two if the observing station is in darkness 
and, in general, observations can be made (subject to 
other checks). 

PROGRAM DESCRIPTION : Data storage is used to input the mean terrestrial 
coordinates (UVW in kilometers) in the array STAOPT and the seleno- 
centric coordinates of the moon (XYZ in kilometers) in the array PM. 
Program formulation then represents the theoretical development given 
in f Papo, 1971], Appropriate checks for observability are indicated with 
comment cards. 

SUBROUTINES REQUIRED : 

A. O. S, U. Project Library: 

1. ROTATE 

2. TRIM 

3. ONEM 

4. REDPI 

B. Fortran Scientific Subroutine Package: 

1. DGMTRA 
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REFERENCES : 

Papo, Haim B. (1971), "Optimal Selenodetic Control, " The 
Ohio State University, Department of Geodetic Science Re port No. 
156. 
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onnooooo 


OPTOBS 


SUBROUTINE OPTDBS « £T tXt JOK, APP^BUNDLE ) 

SUBROUTINE FOR GENERATING OPTICAL OBSERVATION OF TH MCON 

X (1-3) GEOCENTRIC MOON POSITION OR SELENOCENTRIC POSITION OF SATELLITE 
14 - 6) EULEHIAN ANGLES OF MUON SYSTEM 
(7 - 9) EULERIAN ANGLES OF EARTH SYSTEM 
EULERIAN ANGLES : L0^4GITU0E NODE INCLINATION 
X ARE RELATED TO THE ECLIPTIC SYSTEM cv^tcm 

BUNDLE - ANGULAR COORDINATES OF OBSERVATION RaVS IN SYSTEM 


C 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


IMPLICIT ReAL*8CA-H,G-Z) ^ 

01MENSI0NEAB(2) tOBU(3}»TlC3,3)f T2(3t3)»T3{3»3)tSTA0PT(3f3)fPCSP(3) 

DIMENSION SUNC2) t OB (3) f TR(3 I f BV (3) »PCC3) » TM ( 3 13 ) >TB I3t3l t PH '3 #30) t 
PCK(3f30),EUNDLEC2#30) #X(9) 

common /opto/ TB# WTERf SSUN#ETER#PCSP#STI#ST2#ST3#IFLAG 

TH - TRANSFORMATION MATRIX FROM SELENCGRAPHIC TO ECLIPTIC SYSTEM 

PH - SELENOGRAPHIC CARTESIAN COOROINATES OF POINTS ON THE HDON 

PC - SELENOCENTRIC ECLIPTIC CARTESIAN COORDINATES OF PROJECTION CENTER 

TB - TRANSFORMATION MATRIX FROM ECLIPTIC TO B1-B2-&3 SYSTEM * 

OB — UNIT VECTOR OF EARTH OBSERVATORY (GEOCENTRIC ECLIPTIC SYSTEM) 

SUN - UNIT VECTOR OF SUN RAYS IN ECLIPTIC SYSTEM (SUN AT lUFINITY) 

DATA $TA0PT/-*1996,0CSlf ~5042^6961 »3360*7T4a »46B6*1252» ll^6385#4 
P331*0499#5058^2628r 2698 . 0251 ^-27 99 *0019 / 

DATA PM /284.7144 t-164*8233 tl70S.6738 f267»3B28 #742«4968 #134b« 
,298*9995 105O*4336 ,1350^2634 f950*9536 #— 484*8089 #1369* 

P1476 #897*8094 #293.5282 #1457.1998 #297.0151 #1298.3423 #1116. 

P7341 #1013.5746 ,-1080.6066 #894.6916 #1453-2020 #-75.2901 ,946. 

P1413 #1076*7970 #844.2392 #1069.5859 #341.1C61 #-1615.9158 t536. 

P6048 #1506.4731 ,-767.8602 #390.4892 #1577.2501 #544.9516 #470. 

P5874 #341.1861 #1617.9158 #536.6048 #1157.6880 #-1287.5193 #120. 

P7443 #1735.3430 #1*0 ,29.8342 #1106. *♦064 #1322.5432 ,211.3220 # 
PX726.62C7 #-150.2785 #-91.4655 #490.5786 #—1611.7874 #-420.9863 # 
P1499.1838 #733.4201 #-479.5874 #294.7339 #1686.6973 #-302.3192 # 
P1211«2688 #-1091.8823 #-594.9678 #1511.8965 #-131.4927 #-843.1513# 
P1004.3727 #1244.3842 #-679.6325 ,257.4499 #-1221.9682 #-1207.0911# 
P1155.7211 #-614.8371 ,-1140.8012 #1055.5706 #741.8690 #-1163.5211# 
P449.7967 #1120.4735 ,-1250*7900 ,940.5390 «C3.5052 #-1458.1998 # 
P377.2905 #-378.7905 #-1653.5386 #305.4515 #363.0022 #-1671.2764/ 

40K-i”hEANS observatory NO 1 IS OBSERVING AND SAME FOR NOS 2 AND 3 
J0K=9 MEANS THE STATION IS ON A SATELLITE 
00 12 1=1,3 
12 PC(I)=XCI) 


2*P 

I 


14^P 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 


SAN-1.7399359D0+.0172O27913DO^(ET+2498O.D0) 

SSUN=SAN-X ( 4 ) -X ( 5 )-3 . 14 16 
CALL REOP I (SSUN) 

SSUN=$$UN*57.2958 

ETER=SSUN+90. 

HTER=SSUN-90* 
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OPTOBS fConn 


IF{ETER.GT.360.00)ETER=ETER-360.00 
SBN( l)=0C0S(SAN) 

SUNt2)=DSlN(SAN) 

C 

1F«JCK-9J 14,13,11 

14 0015 J=l,3 

15 CBU(J)=STAOPT( J.JOK) 

CALL ROTATE 13, -X f7),Tl) 

CALL ROTATE (1, X t9),T2) 

CALL TRIM tT2,Tl,T3) 

CALL ROTATE (3,- X (E),T1I 
CALL TRIM (T1,T3,T2J 

CALL OMEH (T2,08U,0B) 

DO 10 J=l,3 
10 PC{J)=OBlJ)-PC( J) 

GOTO 19 
C 

13 0016 J=l,3 

16 OBt J »=-PC< J) 

19 CALL ROTATE (3,- X 14), Tl) 

CALL ROTATE (1, X (6),T2) 

CALL TRIM (T2,T1,T3) 

CALL ROTATE (3,- X (5),T1) 

CALL TRIM JT1,T3,TH1 
CALL OGHTRA (TM,T1,3,3) 

CALL ONEM (T1,PC,PCSP) 

C 

EAB(11=DATAN2(-PC(2),-PC{1) ) 

EAB{2)=0ATAN2(-PC(3) ,DSCRT( PC ( 1 ) **2+PC(2 ) ) 

CALL ROTATE 13, EAfi(l),T3) _ 

CALL ROTATE <2,-EAC(2> ,T1 J 
CALL TRIM {T1,T3,T8) 

C 

S8=DS0RT ( OB ( 1 ) ♦♦2+03 ( 2 > **2+03 ( 3 ) ♦♦2 ) 

CHi=<SUNm^afa{n+SUN(2)+03 (2))/SB 
1F{CHI-ST3) 1,1,2 

C CHI - IS IT NIGHT AT EARTH OBSERVATORY 
2 CALL DGHPRO ITH,PM,CH,3,3,30) 

C 

DO 3 1=1,30 
00 5 J=l,3 

5 TR( J)=CH(J,I )-PC( J) 

SC=DSQRT (CMd , I )*=»2+CHl2,l)^^2+CH(3,I )^*2) 

CIJ2={CH(1,I) + SU:U t)+CM(2,I)^SaNI2))/SC 
C CH2 - IS THE MOON POItJT ILLUMINATED 
SU=DSCRT(TRm ♦♦2+TR(2)^*2+TRJ3)^^2) 

CH3=tTR(l )*CM( l,I)+TPv (2)^CM(2,I )+TR(3)*CHC3,I ) )/(SU*SC ) 

C CH3 - IS THf MUON POINT AT LEAST 20 DBG FROM THE LIMB 

CM4=(00( 1)^TR{ l)+crU 2 )^TR (2 ) +LP. (3 )*TR 1 3) )/( SU»SB) 

C CH4 - IS THE MOON POINT AT LEAST 20 DEG ABOVE THE HORIZON OF EARTH OBSERV. 

CALL ONEH (TB,TR,BV) 
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o on o 


OPTOSS fCont) 

R=0SQRT(BVm=»*2^9V(2)+*2+BV(3)**2» 

Bb-NOLE ( 1 1 1 J =DARCCS ( BV U) /R ) 

SUNDLE{2,I)=0ATAN2tBV<3)»aV{2n 

lFtCH2.6T.STl.0R.CH3.GT.ST1.0R.CHA.LT.ST3.0R.0ABS{BV{2)/R )-GT. 
PAPP.DR.0ABSCBV13J/R ) .GT.APP >GOTO A- 
STl - -SINI 5 DEG) SEE SECTION A.61 
ST2 s -SI)4 (20 OEG) 

ST3 = SIN (20 DEG) 

GOTO 3 

4 BUNDLE (It 11=0.00 
BUN0LE(2»I)=0.D0 
3 CONTINUE 
C 

RETURN 
I IFLAG=1 

C IFLAG=1 HEANS THE OBSERVATORY ON EARTH CAN NOT OBSERVE. IT IS DAYTIME 

KR1TE(6,76)J0IC»ET 
RETURN 

U HRITE(6,771ET 
RETURN 

76 FDRHATdHX ///’tSXt* EARTH OBSERVATORY NO 'tlSt* CAN NOT OBSERV 
PE • IT IS DAYTIME SFZO.B) 

77 FDRHAKSX,* IMPROPER VALUE GIVEN TO JQK *iF20.81 
END 
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SUBROUTINE : OUTADJ 


CALL STATEMENT ; OUTADJ (ET, Y, YP, IHL, K, PAR) 

SUBROUTINE PURPOSE : OUTADJ is the program used to adjust the 
numerically integrated physical libration angles to Eckhardt's 
theory TEckhardt, 1970] as described in [Papo, 1971 1 Chapter 3- 

INPUT PARAMETERS : 

A. Scalar (extended precision). 

1. ET is the epoch of interest in Julian days. 

B. Vectors (extended precision), 

1, Y and YP are vectors described in the documentation 
of subroutine FUNPL5. 

C. Variables IHK, K, PAR are dummy variables not used in 
the subroutine, 

OUTPU T PARAMETERS : Variables Parameters are output from the subroutine 
in the common block area MAIOUT in the variables CNORM and CVECT. 

COMMON AREA PARAMETERS : Three common areas MAIOUT, MAIF, and 
ALL are used. The areas MAIF and ALL are described in the 
FUNPL5 Program Description. 

A. Area /MAIOUT/CNORM, CVECT. 

1. CNORM is an extended precision 9 X9 matrix which 
contains the normal matrix formed in the subroutine 
(e.g, each set of observations creates an additional 
’’layer” which increments the normal matrix). 

2. CVECT is an extended precision vector of 9 elements which 
contains the constant vector corresponding to the 

normal matrix in CNORM. 

PROGRAM DESCRIPTION : Section 3.42 fPapo, 1971] describes the theory 
used in developing the subroutine. The subroutine is called by the 
DVDPF series of routine at a desired epoch to evaluate Eckhardt's 
formulation and form the necessary normal matrix and constai.t vector 
for adjustment. Basically the program accomplishes the following 
computations for each epoch of interest. 

A, The differences in the numerically integrated angles and 
angles computed from Eckhardt's theory are computed. 
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B. The differences computed above are transformed as described 
in fPapo, 19711 Equation 3.42. 

C. Appropriate portions of the XJ and Q matrices (Papo's notation) 
are formed from the input Y vector (elements 7 to 42 and 43 to 
60 respectively) as formulated in f Pap o, 19711 Equation 3. 42. 7. 

D. The above computations are then used to increment the normal 
matrix and constant vector (CNORM and CVECT) and control 
is returned to the calling program. 

SUBROUTINES REQUIRED : 

A. 0,S. U. Project Library: 

1. EKHARD’ 

2. ONEM 

B. Fortran Scientific Subroutine Package; 

1. DGMPRD 

2. DGMTRA 


REFERENCES ; 

A. Eckhardt, D. H. (1970). "Lunar Libratidn Tables, " The Moon. 
Vol. 1, No. 2, February. 

B. Papo, HaimB. (1970). "Optimal Selenodetic Control, " The 
Ohio State University, Department of Geodetic Science , 

Report No, 156. 
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OUTADJ 


SUBROUTINE- OUTADJ lETtYtYP,lHLiK,PARl 
IMPLICIT REALMS (A-H,0-2) 

DIMENSION Y(42 ) ,YP ( 42 ) t PAR ( 5 ) »C I o) »PHL{3) f01F{3) tD0(3 J ,EH(3»3) » 
PU(3 >6) tUEM(6,3) iC0W{6,6) , CNORH ( 6» 6 ) , DCW(6 ) fCVECT 1 61 » GOG ( 3 J ,0cN(6J 
PfUUH(3«6}«0UH(6) 

CCHMCN/HAIOUT/CNURH , C VECT ,DEN t RUH/ALL/Q 

DATA EH(3fl).EM(l,l) ,EH(2,IJ,EH(3,3> /I. 00,3*0.00/ 

IFCNUT.E0.l)PAR{5 > = 1.D0 

TEST=<ET-PAR(1) )/PAR(3) 

lF(OMOOfTEST,l. DO). GT. 1.0-05) GOTO 99 

STE*0S1N(Q{3) ) 

SFI=DS1N{0(1) ) 

CFI=OCO$(Q(D) 

EH{1,2)=-SF1*STE 

EHtl,3)=-CFI 

EH{2,2)=-CFI*STE 

EK(2,3)=SFI 

£H(3,2)=DCOS<Q(3) )-l.DQ 
DO I IFL=1,6 
KAP*(1FL-1)*6 
00 1 1=1,3 

1 U(I ,IFL) = (Y(KAP+1 )-Y( 36-^1 ) ) /DEN(IFL) 

CALL EKHAKD ( ET•^244 00a0.D0,0UH,G0G) 

00 2 J=l,3 

2 DIF { J) =Y(36-*-J ) -GOG( J ) *.4848 1368 11095360-05 
CALL CNcH (EH, OIF, DO) 

CALL OGKPRD ( EH,U,UUH ,3 ,3 , 6) 

CALL DGMTRA ( UUH ,UEH , 3,6 ) 

CALL DGKPRO {UEM,UUMtC0W,6,3,6> 

CALL DGKPRO (UEH, DO, DOW, 6, 3,1) 

DO 3 1=1,3 

3 R0H=R0H-^D0 ( I ) **2 
DO 4 1=1,6 

CVECT(1)=CVECT{I)■^0CW(I) 

DO 4 J=l,6 

4 CNORHt I , J )=CNORM (1 , J ) ■>CDW( I ,J ) 

99 CONTINUE 

RETURN 

END 
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SUBROUTINE : OUTEPH 


CALL STATEMENT ; OUTEPH (ET, Y, YP, IHALVE, K, PAR) 

SUBROUTINE PURPOSE : Print routine for output of the intermediate results of the 
simulated ephemeris. 

INPUT PARAMETERS : 

A. Integer dummy variables (single precision). 

1. IHALVE 

2. K 

B. Scalars (extended precision). 

1. ET is the desired epoch in Julian days minus 2440000. 0. 

C. Vectors (extended precision). 

1. Y is an 18-element vector containing the same variables 
as described in the Y vector of FUNEPH. 

2. YP is an 18-element vector containing the same variables 
as the YP vector in FUNEPH, 

3. PAR is a 5-element vector not used in OUTEPH, but 
described in the MAIN program for simulated ephemeris 
generation, 

OUTPUT PARAMETERS : The subroutine prints out ET and selected elements of 
the y vector. The selected elements are the XYZ of the state vector, 
the Eulerian angles of the moon and the Eulerian angles of the earth, 

PROGRAM DESCRIPTION : None 

SUBROUTINES REQUIRED : None 


REFERENCES : None 


Preceding page blank 

■ I 
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o o 


OUTEPH 


SUBROUTINE OUTEPH ( ET ,Y, YP , IHALVE *K ,PAR ) 

WATCH OUT NO K VARIABLE PERMISSIBLE IN THIS SUBROUTINE 
OUTPUT OF SIMULATED EPHEHERIS 
IMPLICITREAL=»=8 (A-HjO-Z) 

DlMENSIONYdS) , YP US ) ,PAR (5 ) 

WRITE(6,70)ET» (YC I) »(Y(I),I=7»9j,tY{I),I=13,15> 

WRITE(l) Y 

RETURN 

70 FORMAT UXt F6. 1 ,3 (1X,3D12 .5 ) ) 

END 
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SUBROUTINE : OUTGAJ 


CALL STATEMENT : OUTGAJ (ET, Y, N2, EM3, Al, AIMS, A3, AZ), 

SUBROUTINE PURPOSE : The subroutine is used to process simulated earth 

based optical observation of the moon for solution of lunar coordinates 
and physical parameters of the moon (initial values of the physical 
libration parameters and .8 , Cgo). The subroutine is called by 
the DVDPF2 routine at epochs where optical bundles have been observed. 

INPUT PARAMETERS : 

A. Scalars (extended precision). 

i. ET is the epoch in Julian days at which the optical 
observation has been taken. 

B. Integer (single precision), 

1. N2 is twice the number of optical rays in the bundle. 

C. Vector (extended precision). 

1. Y is a 60-element vector input to the subroutine by 
either the FUNPL5 or FUNPL6 subroutines. 

D. Matrices (extended precision), 

1. EM3 through AZ are matrices used within OUTGAJ and 
are included in the input parameter set to utilize 
object time dimension facilities of EBM Fortran. 

DATA SET INPUT : The subroutine requires binary input from a computer device 
(either magnetic tape or disk) which in the IBM data device desinator, is 
named unit 3. The input consists of the following parameters: 

A- KOK is the serial number of the bundle of rays to be processed 
(integer, single precision). 

B. ETO is the epoch at which the bundle was "observed** in Julian 
days (extended precision). 

C. 1ST is. the identification number of the bundle. 

D. POS is a 3-element extended precision vector containing the posi- 
tion of the projection center in the simulated selenocentric 
inertially oriented coordinate system (XYZ) in kilometers. 

E. EMBT is an extended precision transformation matrix (3.x 3) 
for coordination transformation from the XYZ to the Bi Ba B3 
system fPapo, 1971] Chapter 2. 

F. ENU and EKA are extended precision 22-eIement vectoi's containiu; 
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the simulated optical observations. 

G. IPO is the identification number of the control point 
being observed. 

DATA SET OUTPUT : As described below, the subroutine essentially forms 
and increments a normal matrix and constant vector for adjustment 
of optical observations. The "output" of the routine (e. g. the incre- 
mented normal matrix and constant column) are output on a data set 
device identified as K3 which is provided through the common area 
MAIOUT. The output is written in binary mode. 

COMMON PARAMETERS : Common areas are used extensively to transfer data 
in and out of the subroutine. 

A. Area /MAIOUT/POS, EMBT, ENU, EKA, IPO, TJ, SIG6, 

SIG7, SGGZ, K3, K2, 1ST, MM. 

1. The parameters POS through IPO are described in the 
Data Set Input section. 

2. TJ is an extended precision matrix of (3 X 22) elements, 
and contains the Cartesian coordinates of the control 
points in the selenodetic coordinate system (in km). 

3. SIG6 is an extended precision 3X3 array containing 
the a priori covariance matrix of the orientation 
elements of the optical bundle of rays (in radians 
squared). 

4. SIG7 is an extended precision scalar which is the 

a priori variance (in radians squared) of the optical 
observations relative to the Bg Bis coordinate system. 

5. SGGZ is a 3 X 3 extended precision array containing the 
covariance matrix of the projection center (km squared) 
in the selenocentric inertially oriented system. 

6. K3 is a single precision integer which identifies the 
number of the data set device to be used for output. 

7. K2 is a single precision integer identifying a data set 
device to be used for temporary storage. 

8. 1ST is a single precision integer which is not used in the 
program. 

9. MM is a single precision integer containing the number 
of the first optical bundle to be processed from bundles 
sequentially stored on unit 3. 

B. Area /dVOUT/ETO, K3\l, ICW, N, ISKIP. 

1. ETO is an extended precision variable containing 

the epoch in Julian days of an optical bundle read from 
unit 3. 
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2. KM is a single precision integer containing the total 
number of bundles stored on unit 3. 

3. ICW is a single precision integer containing the 
sequential number of the optical bundle being processed. 

4. N is a single precision integer containing the number 
of rays in the next bundle. 

5. ISKIP is a sibgle precision integer used as a flag to 
indicate that some of the sequentially stored optical data 
are not to be processed rPapo, 1971] Page 204. 

C. Area /ALL/. The common area /ALL/ is input from the subrouc:.- 
FUNPL6 and the parameters are described in FUNPL5 documen- 
tation. 

PROGRAM DESCRIPTION : The subroutine formulation follows the theoretical 
development given in f Papo, 1971], particularly Sections 2, 52 and 2. 6 
(Pages, 59, 68-70). Basicall}' the logic of the program can be divided.m: : 
four- steps; 

A. The calling subroutine DVDPF2 supplies the current physical 
libration angles and the state transition and parameter sensitiviv ' 
matrices as input parameters. 

B. From a previous step of integration (e.g. the last call of OUTG.’- 
the optical data and auxiliary information has been read from 
unit 3 and stored in common storage. 

C. Calculations according to the theory of fPapo, 1971] Chapter 2 
are accomplished and the increment of the normal natrix, etc. 
are added to the accumulated results and stored on a disk unit 
(direct access device). 

D. The optical and auxiliary data of the next, sequentially stored, 
bundle of rays is read from unit 3 and control is returned 

to the calling subroutine DVDPF2. 

SUBROUTINES REQUIRED : 

A. O. S.U, Project Library: 

1. REDPI 

2. ROTATE 

3. TRIM 

4. ’ LUCA 

5. ONEM 

B. Fortran Scientific Subroutine Package: 

1. DGMTRA 

2. DGMPRD 

3. DMINV 


REFERENCES : 

Papo, Haim B, 0971). "Optimal Selcnode He Control, " 
The Ohio State University, Department of Geodetic Science , 
Report No. 156. 
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OUTGAJ 


OUTGAJ 


SUBROUTINE OUTGAJ {ET,V ,N2,£H3iAl»AIH3,A3,AZ) 

1KPLICITRSAL*8 lA-H,0-2) 

DIHENSIOM Y{60J,SEV(6),2(3),UQ«3,9),T£1(3,3J,TE2{3,31,TE3C3,3) 
PfTE«l(3*3) tEHHi3f3> ,EHB (3 »3) ,£M3T(3 1 3 ) , EH3{.\'2»N2) yAl (N2.75) , U ( 3) 

P, rPQ(2i)tENUJ22J»BKAC22)TU5(3},TJ{3,22),T£PAr3t3),P0S(3),EHJH{3f3J 
P, EHT{3lf86(2,3)tETT(2,3>,T£M2l2f3}tS23t2f9),TEH5{2,31,B10(2t3) 

P» SlGZ(3T3)yEHZ(2 «2 ) (2 »3 } tW < 2 ) t A 3(N2, 3),AZ(N2« 2} ,LUH( ,aOL 
PJ+9I »SlG6{3f3>tAlH3{75,N2),UU(A4),TrH3(3),VEN2{44),VEC(75),Nt31) 
CtMMON/MAICUT/ POS,E«BT,£NU,EKAfTJ,SlG6tSIG7»SlGZ,K3tK2TlPOf ISTfHH 
P/DVOUT/ ETO»KH ,ICW,N,1SI<IP /ALL/ SEV,RO«»NUT 
DATA PI /3. 14159265358979300/ 

. FIICW.LT.HH) GOTO 17 
HRITE(6,73) ET,ETO, ICW, N1 ,N2 
IF(ICH.GT,K.M1. GOTO 99 

C BEGIN CALCULATIONS FUR THE PRESENT BUNDLE OF DIRECTIONS 
00 1 J=l»9 
JJ=J#6 
00 2 1=1,3 
2 UQ{I,J)=YtJJ+n 
1 UQtl,J)=UQ{l,J)-UCC2fJ) 

C UO IS THE MATRIX OF PARTIALSFOR PH. L. tTRANSITION AND SENSITIVITY) 
2tl)=Y<l)-Y(21+SEV(l)-SEV(2)+Pl-HSEV(3)-SEVC4))*(bT-222.5Li)) 

CALL REOPl (Ztin 

2t2)=Y(2>+SEV(2)»SEV(A)*EET-222.500) 

CALL REOPl (Z(2)) 

Z(3)=Y(3)+.026769D0 

C 2 ARE THE EULERIAN ANCLES OF THE MOON 1 - PHI 2 - PSI 3 - TETA 

CALL ROTATE C3,-Zm,l£l) 

CALL ROTATE ll,Z(3),TE2) 

CALL ROTATE 13,-212) ,TE3) 

CALL trim tTE2,TEl,TLHl) 

CALL TRIM (TE3,TEM1,EHHJ 

C EHM IS THE TRANSFORMATION MATRIX FROM SELENOOETIC TO INERTIAL 

CALL DGKTRA I EPBT tEHB ,3 ,3 ) 

CALL LUCA (l,-l.OC,TEl) 

CALL TRIM tTE3,TEl,TE2) 

CALL TRIM ITE2,TEM1,TE3) 

CALL LUCA 13,1.00,TEHI) 

CALL TRIM ( EHM,TE.'U,Tcl) 

CALL TRIM <TcHl,EHH,TE2) 

C TEltT£2,T€3, ARE NEEDED FOR THE MATRIX OF PARTUALS B23 
DO 12 1=1,66 
DO 12 J=1,N2 
12 AHJ,I)=0.00 

INITIALIZING PART OF THE A1 MATRIX 

THIS COMPLETES MATRICES COMMON TO ALL THE POINTS IN THE BUNDLE 
L=N2/2 
DO 13 1=1, L 
H=1P011) 

SN=DSINUNu(I) ) 

CN=OCOS(ENU(D) 
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n o 


OUTGAJ (Cont) 


SK=OSlNI£KA(in 

CK=DC0S<6KA(IJ ) 

U{1)=CN 
U(2)=SN^CK 
UC3)=SN*SK 

C U ARE THE DIRECTION COSINES OF RAY (I) IN THE B1B2B3 SYSTEM 
CALL ONEH JEHB,U,UB) 

00 3 11-1,3 

TEH4{II,l)=TEini .D^TJdfMJ+TEKlI.Zl+TJtB.HJ+TEKIltSI^TJOtH) 
TEH4{II,2)=TE2(II ,l)*TJ{i,M)+TC2ai,2)*TJ{2,M)+TE2ai,3)»TJ(3,H) 
TEM4(II,3)=TE3tII ,l)*TJ(l,M) + T£3{ri,2)*TJ(2,H)+TE3(II,3)=i'TJ{3,M) 

TCH4 IS THE MATRIX IN BRACKETS FOR B23 
TJ ARE THE "08SERVCO” COORDINATtS OF THE HOON POINTS 
TEH3(11)=-TJC1I,/I) 

DO 3 J=l,3 

EHUHCII,JT=UBfII)*U0( J) 

I F( 1 1, £0. J ) EMUH ( II , 1 1 ) =EHUH ( II , II J-1.00 
3 CONTINUE 

C EHUH IS THE EXPRESSION IN BRACKETS FOR Bl { NEGATIVE) 

CALL ONCH IEHH,TEM3,EHT) 

DO 11 J=l,3 

11 EHT{ J)=EHT(J)+PUS (J) ' 

C EMT IS POSITION OF PROJ. CENIR WITH RESP. TO OBSERVED PNT IN INERTIAL COMP 

RO=OSCRT(EHTtl>**2+tMT<2)=«'=«'2+EMT{3)**2J ' 

B6( l,l)=-CN»(CK*Uf3)-SK*U(2) ) 

B6( 1,2)=-SN*U13)-CN*SK-U( 1) 

86(1,3)= SN*U(2)+CN»CK*U(1) 

66(2,1)= SK*U(3)+CK«‘U(2) 

B6{2,2>=“CK=»Utl) 

B6(2,3>=-SK*U(1) 

C B6 IS THE MATRIX OF PARTIALS FOR THE E ORIENTATION PARAMETERS 

ETT( 1,1)=-SN/R0 
ETT(1,2)= CN-»CK/RO 
ETT(1,3)= CN*SK/R0 
ETT(2,1)= 0.00 
ETT(2,2)=-SK/R0 
ETT(2,3)= CK/Rt> 

C ETT IS THE ET/RO MATRIX 

CALL DGHPRD ( ETT, EMBT ,T£M2, 2 ,3 , 3 ) 

C TEM2 IS ET^Ht)T/RO 

CALL DGMPRO ( U M2 , TEH4, TE H5 , 2,3,3 ) 

CALL DGHPRD (Tl M5 ,U0, B23, 2,3 ,9 ) 

CALL DGMPRD ( TEM2 ,EHUM,6 10, 2 ,3 ,3 ) 

CALL DGHPRD ( B 10 , SIG2 ,T£H5, 2,3, 3 ) 

DO 4 K=l,2 
DO 4 J=l,2 

A EHZ(K, J)=TEM5(K ,1 )*B10( J , 1) +TEM5 ( K,2 )»8 10 ( J ,2 )+TEM5 (K,3)*S10 ( J, 3) 

CALL DGMPRD ( 610, EMM, 81 , 2 ,3 ,3) 

C 81 IS THE MATRIX OF PARTIALS FOR THE TRIG POINTS ON THE HOON 
DU 6 K=l,2 
H(K)=0,D0 
DO 5 J=l,3 
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OUTGAJ (CoDt> 


BKK t J>=BI(K,J)*20626A.8062DO 
&10{KtJ}=810(KrJ)*206264. 806200 

5 HCK)=H(K)+TEM2(K,JS«(UB« J)+fcHT(J)/RQ> 

6 W(K)=H(KJ*R0»206266.8C6200 
H IS TH5 FREE TERM 

CALCULATIONS CF 61 » B23 t B6 t BIO ARE COMPLETED 'FOR POINT H 
THE SUBMATRICES CREATED ARE STORED IN A1 , A3 , UU AND A2 

DO 10 K=1t2 
1I=(I-1)»2*K 

uuni)=w(K) 

OD 7 Jsl,3 
A3ClI,J)=B6(KtJ) 

7 A1(II,JJ+J1=B1(K,JJ 
DU a J=l,9 

8 Al(II,66-»J)=B23CKtJ) 

DO 9 J=l,2 • 

9 A2m,J)=EMI(K,J) 

10 CONTINUE 

13 CONTir4UE 
HRITE(6,74> UU 
DO 15 1=1, N2 
DO 15 J=1,N2 
EH3n,J)=0,00 
00 14 K=l,3 

14 EH3(I,J) = EH3n,J)+A3(l,K)*(A3(JTl)*SIG6CK,ll+A3(J,2)*SIG6lK,2) + 
PA3(J,3)*SIG6(K,3} ) 

IFU.EQ.J) EM3(1,J)=EH3(1,J)+SIG7 

15 CONTINUE 

HR1TE(6,74HEH3<I,I),I=1,N2) 

FAT=O.DO 
DO 16 L0=1,N2 

16 FAT=FAT+DLQG10(Eii3(LQ,L0) ) 

FAT=1.D0/10.00**(FAT/0FL0AT«N2) ) 

DO 22 L0=1,N2 
00 22 LP=1,N2 

22 EH3(LO,LP)=EH3(LO,LP)*FAT 

CALL DHINV {EM3,N2,DE7,LOM,MOL> 

DO 24 L0=1,N2 
DO 24 LP=1,N2 

24 EH3(LO,LP»=EM3(LQ,LP)*FAT 
TEST=O.DO 
DO 25 L0=2,N2 
LC=L0-1 
DU 25 LP=1,LC 

25 TBST=TESH-(tH3(LP,LO)-EH'3CLO,LP) )**2 
KRITE(6,74) DET, FAT, TEST 
WR1TE(6,74) (EK3<I,I J ,I=1,N2) 

C MATRIX M3 IS INVERTED 
DO 19 1=1,75 
DO 19 J=1,N2 
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OUTGAJ (Cent) 


AlM3tItJ)=0.D0 
DO 18 Ksl,N2 

AlM3{I,J)sAiM3(I,J)+Al(K,I)»EM3(K,J) 

CONTINUE 
REWIND K3 
REWIND K2 
DO 21 1=1,75 
READ(K3J VEC 
DO 20 J=l,75 
00 20 K=1,N2 

VEC(J)=VEC{J)+A1H3(I,K}’*A1<K,J) 

KRITEtN2) VCC 
READ (K3) VEC 
00 23 J=l,75 
00 23 K=I,N2 

VEC( J)=VEC( J}+A1H3( J,K)»UU(K) 

HRITE(K2) VEC 
DO 28 1IK=1,75 
READ(K3) VEC 
DO 27 H=1,L 
KA1={K-1)*2+1 
KA2=(H-TJ*2+2 

VEN2tKAl)=AlH3ClIK,KAl)*A2{KAl,l )+AlH3(IIK,KA2)*AZIKA2,l) 

VEN2 (KA2 J=A1H3(IIK,KA1 l*A2(KAl,2)+AlH3(IIK,KA2)t:A2{KA2,2) 

CONTINUE 
DO 29 Hl=l,75 
00 29 H=1,N2 

VEC(Ml)=VtC(Ml)+VEN2(*H>*AlH3{Hl,MJ 

WRITF(K2) VEC 

CONTINUE 

IAX=K3 

K3=K2 

K2=1AX 

IFnCH-GE.KH) GOTO 99 
IFtlSKIP.EQ.O) GOTO 17 
DO 707 1=1,1SK1P 
REA0(3,EN0=99) KOK 
ICW=ICW+1 
Nl=NtICW+l) 

READ(3,END=99) KOK, £TQ, I ST,P0S,EM9T, flPOI I J ,ENU( II ,EKA (I ) , 1=1 ,Ni) 
HRITE(5,73J ET,ETO, ICK.Nl ,N2,K0K 
lFnKOK-lCH).NE.l) GOTO 98 
RETURN 

WRITe(6,7lJ ICW.KDK 

NUT=1 

RETURN 

NUT=l 

RETURN 

F0RMAT(//10X,2I5) 
format (/1>X,2F12. 6,415 ) 

FORMAT (/<5X,12U9.2) ) 

END 
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SUBROUTINE : PL ADIS 


CALL STATEMENT : PL ADIS (ET, EX, XM, XCE, X, D, PE, PM, EP, 

EL, SE, SM, PSM, Bl, B2) 

SUBROUTINE PURPOSE : Computation of partial derivatives of simulated 

laser distances with respect to the considered parameters as outlined 
in r Fajemirokun, 1971] Chapter 5. 

INPUT PARAMETERS : 

A. Parameters ET, XE. XM, XCE, X, D, PE and PM are 
described in the subroutine LADIS Program Description. 

B. EP is a 6-element extended precision vector representing 
the integrated Eulerian earth angles 

fFajemirokun, 1971] Page 164. 

C. EL is 6-element extended precision vector representing 
the integrated Eulerian moon angles 

rpajemirokun, 1971] Page 192, 

D. SE is a 3x6 extended precision array which contains the 
state transition matrix of the earth system as described in 
fFajemirokun, 1971] Page 64, 

E. SM is a 3x6 extended precision array containing the state 
transition matrix for the moon. 

F. PSM is a 3x3 extended precision array containing the state 
transition matrix of the dynamic parameters of the moon, 

OUTPUT PARAMETERS : 

A. Bl is a Z1 element extended vector containing: 

1. Partials with respect to the earth station position (1-3). 

2. Partials with respect to the earth Eulerian angles (4-9). 

3. Partials with respect .to the lunar station position (10-12) 

4. - Partials with respect to the lunar Eulerian angles (13-lS). 

5. Partials with respect to the dynamic (Csot ^) 
parameters of the moon (19-21). 

B. B2 is a 3-element extended precision vector containing the 
partials with respect to the geocentric coordinates of the 
center of mass of the moon. 

PROGRAM DESCRIPTION : The program computes the partial derivatives 
formulated in ‘^Fajemirokun, 1971] and in the order giv'en above. 

Preceding page Uank 
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^SUBROUTINES REQUIRED : 

A. O. S. U. Project Library: 

1. LUCA 

2. RATATE 

3. ONEM 

B. O. S,U. Utility Library; 

1. MMULT 

2. MSCALE 

REFERENCES ; 

Fajemirokun, F. A. (1971). ’’Application of New Observational 
Systems for Selenodetic Control, ” The Ohio State University, 
Department of Geodetic Science, Report No. 157 . 
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PIADIS 


SUBROUTINE PLAOISI ET , X£ »XH, XCE, X,D,PE ,PM, tPf EL« SE. SM, PSM, Bit 62) 
lMPLlCirR6AL=!=a( A-H,0-Z) 

DIMENSION X£{3) tXM(3 » ,XCt (3 J tX(3) tPb(3t3)»PM(3t3 ) tEP( 6J »£L16) , 
»&2(3J,T1(J) ,SEC3,6) ,SH(3,6),PSH(3,3J,TEH1I3.,3J tTEl(3t3)* 
♦TE2t3t3) ,TE3ti.3) ,T/t3) ,TfcH2 (3 ,5) t 6 C 6 ) 

OrHENSION A(3) tP(3)tCl(3t3i .02(3.3) 

DR=1.00/D 
D020I-1.3 
0020J=i.3 
C1H,J)=0.D0 
C2U.J)=0.00 
20 CCN7I;<.UE 

C1(1,1)=1.D0 
Clt2.2)=1.00 
C113.3>=1.D0 
C1(3.2)=-1.D0 
C2(1.3>=1.D0 
02(2.21=1-00 
C2{3, i)=l-CiO 
02(3.2)=-!. 00 

C START FCRMINO PARUaLS TJ.R.T.THE GEOCENTRIC COORDS. OF EARTH STN. 

0ALLMMULT(X,Pfc,l,3,3. U) 

CAI.LMSCALE(-DR,T1.1,3) 

00301=1,3 
Dl(n = TllI) 

30 CONTINUE 

C THIS COHPLtTtS PARTlALS W.R.T.GEOC .COOPOS.'OF EARTH STATION 

C START FORMING PARTIALS ’ri.R.T. EARTH EULtRlAN ANGLESt INITIAL CONO-) 

0 ALLLUC A ( 1 , 1 .CO , TEM 1 ) 

OALLROTArE{3,-bP(2» , TE2) 

CALLRQTATEU.FPID.TEl) 

CALLR0TATE(3,-EP(3) ,T£3) 

CALLCiMEHiTcj.XL.Tl) 

CALL(JNtH(TEi.Tl, T2) 

CALLQNEH( TEMl.T2.il) 

CALL0NEH(TE2,T1.T2) 

A(1)=0.DC 

00331=1,3 

A(1)=A(1)+X(I)*T2(I) 

33 CONTINUE 

A(l)=AU)/0 

CALLLUCA(3.1.D0,TEH2) 

CALLONtMCPE.XE, Tl) 

CALL0NEH(TEH2,T1,T2) 

A(2)=0.D0 

00361=1,3 

A(2)=At2)+XtI )*T2( I ) 

3o CONTINUE 

A(2)=-A(2)/D 

CALL0NEH(TE.'i2,X£,T1) 

CALLONEKIPE.ri, 12) 
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PLADB fContt 


A(3)s0.U0 

0D39I=l,3 

A(3)®A{3»*X<1)»T2U) 

39 CONTINUE 
Al3»a-A*3I/D 
CALLHKULTUtClf l>3»3fTI 
CALLHHULTIF,SE,l*3,o,B) 

D040I»i#6 

40 CONTINUE 

C THIS COWLETES PARTlAUS W.R.T.INITIAC CONOS. FOR EARTH»S ORIENT. 

C START FORKING FARTlALS FOR S6LEN0O. COORDS, OF LUNAR POINT 

CALLHMULrtX.PM,l,5,3,TlI 
CALLMSCALE(DR,TI,1,3) 

D0501»l,3 
BUI^9»=T1III 
50 CONTINUE 

ENDS FORMATION OF PARTIALS H.R.T.SELENOO. COOROS. OF LUNAR POINT 
START FCRMli»& PARTIALS W-R.T. HQON'S EULERIAN ANGLS.UNIT. CC’NDS.J 
CALLRCTATEU,£L(1> .TEll 
CALLRUTATtI3t-EL12),Te2) 

CALLROTATF{3,-eL(3J,Te3i 
CALLONtMI TE3,XM,T1) 

CALL0N£M{T61,Ti,T2) 

CALL0NEHtT£Hl,T2TTl» 

CALL0NFHlTE2tTl,T2» 

ACH=O.DO 
00631=1,3 

Am=Am+XIII4T2UI 
63 CONTINUE 
AI1)=A(1J/D 
CALLONeMlPM,XH,Tl> 

GALLON EM 1TEH2, n.T2) 

A{2)=0.00 
00661=1,3 

A(2)=A(2|i»^X|U*T2(l| 

66 CONTINUE ' 

A<2)=-A(2l/0 
CALLQNEK(TEH2,XH,T1> 

CALLCjNEM(PH,TI,T2} 

At3)=0.L0 
00691=1,3 

A13I»AI3I+X{1I*T2(1> 

69 CQNllNUE 
A(3I=-AI3)/0 
CALLMHULT|A,C2,1,3,3,F1 
CALLKKUH I F,SK, 1,3,6,61 
00701=1,6 
Bin+12)=B<II 

70 CONTINUE 
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PLADIS (Cont) 


THIS ENOS FORMATION OF PARTIALS OF MOON’S ANGLES UNITIAL CONOS.) 
START FORMING PARTIALS FOR DYNAMIC PARAMETERS OF MOON 
CALLHMULnF,PSM,l,3,3»bJ 
ooaoi=i,3 

80 CONTINUE 

ENDS PARTIALS FOR DYNAMIC PARAK. OF HUoN — C22,C20,E£TA. 

FORMING PARTIALS FOR CCCROS (GEOCENTRIC! £jF‘ LUNAR CENTER 
00901=1,3 
S2(n=X(I)/0 
90 CONTINUE 

ENDS PARTIALS FOR COORDS OF LUNAR CENTER 
ENDS FORMATION OF ALL PARTIALS 

VALUES RETURNED IN B1 AND 82 MATRICES WILL HAVE TO BE PUT 
IN PROPER LOCATIONS IN MAIN PRQGRAM/OUI PUT SUBROUTINE. 

999 RETURN 
END 
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SUBROUTINE : PM AT 


CALL STATEMENT : PM AT (PSI, THETA, PHI, P) 

• SUBROUTINE PURPOSE : The subroutine computes the transformation matrix 
necessary to rotate an earth fixed or lunar fixed coordinate system to 
. the celestial system. 

INPUT PARAMETERS : 

A. Scalars (extended precision). 

1. PSI is the first Eulerian angle i!) in radians. 

2. THETA is the second Eulerian angle 0 in radians. 

3. PHI is the third Eulerian angle <0 in radians. 

OUTPUT PARAMETERS : 

A. Vector (extended precision). 

1. ' P is a 3 X3 array which contains the transformation 
matrix to rotate from the fixed body system to the 
simulated earth-moon system. 

PROGRAM DESCRIPTION : The output array (P) is formed by multiplying 
successive rotation matrices : 


P = Ri(e) R3(-<P) 

SUBROUTINES REQUIRED : 

O. S. U. Project Library: 

1. ROTATE 

2. TRIM 

REFERENCES : 

A. Fajemirokun, F. (1971). ’^Application of New Observational 
Systems for Selenodetic Control, ” The Ohio State University, 
Department of Geodetic Science. Report No. 157. 

B. Papo, Haim B. (1971). ’’Optimal Selenodetic Control. ’’ 

The Ohio State University, Department of Geodetic Science , 
Report No. 156. 


Preceding page blank 
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SUBROUTINE PMAT (PSI, THETA, PHI, P) 

THIS ROUTINE COMPUTES THE P-MaTKIX NECCESSARY TO ROTATE FROM AN 
EARTH-FIXED OR LUNAR FIXED COORD. SYSTEM TO AN “'INERTIAL** SYSTEM 
INPUT QUANTITIES ARE EULER ANULES PSI , THETA , AND PHI 
P=R3(-PSI) .RH THETA ).R3( -PHI} 

IMPLICIT REAL=^8<A-H,Q-Z) 

DIMENSION P(3,3) ,Ti(3,3) ,T2{3,3I,T3(3,3} ,TA(3,3) 
CALLR0TATE(3,-PSI,T1) 

C ALLROTATE (1, THETA, T2) 

CALLR0TATEC3,-PHI,T3) 

CALLTftlH{Tl,T2,T4) 

CALLTRIH(T4,T3,P) 

RETURN 

END 
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SUBROUTINE : POLE 


CALL STATEMENT : POLE (DT. XPM, YPM) 

SUBROUTINE PURPOSE : Computation of the coordinates of the true pole from 
the pole of the Conventional International Origin using values from the 
International Polar Motion Service (IPMS) from 1958 to 1970. 5. 

INPUT PARAMETERS : 

A. Scalar (extended precision). 

1. DT is the epoch of observation in Julian days 
minus 2400000. 5 days. 

OUTPUT PARAMETERS : 

A. Scalars (extended precision). 

1. - XM is the polar motion coordinate in the x direction 

as defined b}" the IPMS in seconds of arc. 

2. YM is the polar motion coordinate in the y direction 
as defined by the IPMS in seconds of arc. 

PROGRAM DESCRIPTION : The program uses polar motion values stored in 
data statement storage. The program logic is as follows; 

A. The input date is checked to see if the date falls within the 
table limits (1958-1970. 5). If not, an error message is printed 
and control returned to the calling program. 

B. If the input date is valid, a linear interpolation of the tabulated 
polar motion values is accomplished. 

SUBROUTINES REQUIRED ; None 

■ REFERENCES ; 

Fajemirokun, F. A. (1971). "Application of New Observational 
Sj'stems for Selenodetic Control, " The Ohio State University, 
Department of Geodetic Science, Report No. 157. 


- 115 - 


POLE 


SUBROUTINE POLE(OT, XPM»YPH) 

DOUBLE PRECISION DT,XPM»YPM 
DIMENSION PH(251,2),PKT(2) 

DIMENSION PMXK90) »PHY1{90) ,PMX2< 90 ) »PMY2 1 90 ) ,PMX3 1 15 J tPMY3( 15) 
EQUIVALENCE! PM( If 1) iPHXUl) ) » (PH( it2 ) fPMYH 1) ) , (PM C9l»l) f PHXZt 1 ) ) t 
1 ( PMC 91,2) fPHY2 (1 ) ) , (P«{181,l),PHX3(l) ) TCPHnSliZhfPHYSd)) 
DIMENSION PMX4C5) ,PMYA(5) 

equivalence (PH( I9e,,i) ,PMX4U) ) , (PM(196,2) ,PHY4(1) ) 

DIMENSION PHX5(51] »PHY5(51) 

EQUIVALENCE (PH(201, 1 ) ,PHX5 { 1) ) , (PHI 201 ,2 ) ,PHY5(i I ) 

POLAR MOTION TABLES FURNISHED BY lOHLINSOH (TAKEN FROM IPMS) 

DATA PMXl/ 


2 

3 


0 . 111 , 

0.218, 


0.188, 0.237, 0.348, 0.398, 0.398, 


5 0.308, 0.296, 0.261, 0.202, 0.135, 0.073, 

6- 0.026,-0.072,-0.096,-0.107,-0.103,-0.087,- 

7 0.070, 0.080, 0.109, 0.117, 0-117, 0.109, 

8 0.064, 0.062, 0.057, 0.046, 0.034, 0.030, 

9 0.042, 0.041, 0.059, 0.028, 0.019,-0.010,- 
A 0.008, 0.027, 0.C47, 0.071, 0.095, 0.120, 

DATA PHX2/ 

1 0.171, 0.157, 0.128, 0.096, 0.C56, 0.C17,- 
2-0.110,-0.121 ,-0.119,-0.105,-0.076,-0.038, 
3 0.191, 0.239, 0.274, 0.301, 0.281, 0.237, 

4- 0. Oil, -0.069, -0. 122, -0. 171, -0. 206, -0.196, • 

5- 0.055, 0.006, 0.076, 0.166, 0.216, C.260, 

6 0.250, 0.219, 0.161, C-099, 0.062,-0-012, 

7- 0. 165,-0.196,-0.196,-0-176,-6.130,-0.072, 

8 0.168, 0.201, 0.221, 0.227, 0.220, 0.196, 

9 0.000,-0.029,-0.058,-0.086,-0.105,-0.116, 
A-0. 066, -0.057,-0. 010 , 0.052, 0.096, 0.117, 

DATA PHYl/ 


0.606, 0.610, 
0.063, 0.007,- 
0.182, 0.209, 
C.271, 0.269, 
0.160, 0.155, 
0 . 212 , 
0.271, 

0.083, 
0.376, 
0.091, 
0.230, 
0.339, 


1 0.022, 0.C98, 

0.187, 

0.265, 

2 0.676, 0.667, 

0-611, 

0.365, 

3-0. 007, -‘0.038,-0.05 7,- 

-0.0tj6,— 

6 0.285, 0.360, 

0.372, 

0.393, 

5 0.260, 0.201, 

0^X43t 

0.090, 

6 0.059, 0.096, 

0.123t 

0.153, 

7 0.3COr 0.306, 

0.301, 

0.288, 

8 0.150, 0.151, 

0.158, 

0.161, 

9 0.156, 0.157, 

U. 165 , 

0.176, 

A 0.309, 0.316, 

0.312, 

0.306, 

DATA PHY2/ 

1 0.132, 0.092, 

0.068, 

0.060, 

2 0.200, 0.268, 

0.295, 

0.329, 

3 0.369, 0.307, 

0.251, 

0.193, 

4 0.005, 0.041, 

0.078, 

0.120, 

S 0.655, 0.667, 

0.459, 

0.626, 

6 0.123, 0.085, 

0.060, 

0.066, 


0.191, 

0.290, 

0* C67 f 
0.356, 
0.139, 
0.163, 


-0.097,-0.032, 

0.0J6, 

0.368, 0.350, 

0*280, 

-O . 160 » — O . 1 5 3 » - 

0.151, 

0.223, 0.272, 

0.299, 

C.066, 0.035, 

0.U12, 

-0.039, 0.006, 

0.040, 

0.092, 0.076, 

0.065, 

0.032, 0.060, 

0.043, 

-0. 027,-0. 021,- 

-0.009, 

0.166, 0.162, 

0.173 

-0.019,-0.056,- 

-0.086 f 

0.009, 0.070, 

0.134, 

0.176, 0.112, 

0.C48, 

-O, 169t-K»«139f^0-iClt 

0.261, 0.239, 

0.255, 

-0.0o7,-G.120,- 

-O.lbO, 

-0.003, 0.071, 

0.127, 

0.138, 0.075, 

3.033, 

-0.119,-0.115,-0.106, 

0.125, 0.123, 

0.115/ 

0.663, 0.678, 

0.493, 

0.165, 0.097, 

0.C43, 

0.032, 0.120, 

0.211, 

0.601, O.370, 

0.320, 

-0.012,-0.007, 

0.G25, 

0.238, 0.263, 

0.288, 

C.220, 0.189, 

0.161, 

0.153, 0.150, 

0.151, 

0.262, 0.276, 

0.297, 

0.246, 0.214, 

0.175/ 

0*104, 0.128, 

0.16U, 

0.383, 0.337, 

0.375, 

0.046, 0.008, 

-0.020, 

0.294, 0.353, 

C.^X2, 

0.275, 0.219, 

0.168, 

0.069, 0.103, 

0.153, 


COCO 

0001 

0002 

0003 

0004 
0003 
0006 
0007 
oooa 


CO 10 
0011 
0012 
CO 13 
CO 14 
0013 
0016 

0017 

0018 

0019 

0020 
0021 
0022 
0023 

0025 

0026 
C027 
0028 
0029 

C030 

0031 

0032 
00^3 
003 ^ 
0035 
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7 

0.226, 

0.286, 

0.334, 

0.374, 

0.408, 

0.434, 

0.444, 

0.433, 

0.399, 

CO 5 6 

8 

0.349, 

0.303, 

0.259 1 

0.221, 

0.186, 

0.156, 

0.131, 

0.114, 

0.103, 

0037 

9 

0.090, 

0.100, 

0.108, 

0.124, 

0.149, 

0.161, 

0.215, 

0.255, 

0.29b, 

0038 

A 

0.330, 

0.344, 

0.345, 

0.337, 

0.324, 

0.308,' 

0.291, 

0.273, 

0.253/ 

0039 

DATA PMX3/ 









1 

0.099, 

0.079, 

0.056, 

0,031, 

0.012, 

-0.001, 

-C.006, 

— O.C08 , 

-0.002, 


2 

0.0X2, 

0.035, 

C.055, 

0.046, 

0.027, 

0.008 

/ 





DATA PMX4 /-'0*010,-0,029,-0*049»~0*063,-C.066/ 

DATA PHY3/ 

1 0*233t 0.213t 0.194t 0*177, 0.165 t 0.157, 0^155, 0.154, 0.152, 

2 0.156, 0.163, 0.172, 0.163, 0.195, 0-208 / 

DATA PHY4 / 0.220, 0.234, 0.249, 0.269, 0.289/ 

DATA PHX5/ 

1“. 056, -.037, -.014,. 008, .031 , .05 1 , .064, .067, . 064, .060 , .088 , . 119, .1 0 
26 ,.054,. 008, -.027, -.056, -.034,-. 109,-. 123 127, 120 102 

3, -.073 ,-.033,. 010,. 052,. 091,. 125,. 154,. 174, .185,. 164,. 168, -127,. 07 
47 , .02 9, -.021, -.071, -.115, -.157, -.184,-. 184, -.166, -.135, -.100, -.06 
53 ,-.025, .017, .083, .154/ 

DATA PMY5/ 

1.302.. 308.. 308.. 302. .290 . .276 . .260. .245 . . 23 1 . .2 16. .202 . . 183. . 166, 

2. 157. . 156.. 161.. 172, . 197 , . 233, .265 , . 2R9 , .310, .330 , .350, . 370 , 

3. 386. . 392.. 386.. 367. .337.. 502.. 260.. 212.. 167.. 134.. 115.. 105. .104, 
4.114,-134, .168,. 2 16, .273,. 333,. 384, .419, .449, .465 , .463 , . 436 , .39 1/ 

A=:(OT-0. 36203861135)^0. 54758 185D-1 
L=A + 1.0 

IF(L.LT.2)G0T0901 

IF(L.GT.198)G0T0901 

TL==t 

AN=A+1.C-TL 

B=AN + (AN-1.0)/4‘.0 

D010I=1,2 

0£L0=PH C L , I ) -P H < U-1 , I ) 

0EL1=PH(L+1,I)-PH(L,I) 

DEL2=PM ( L+2, I ) -PH (t+ 1,1) 

PMT(I)=PH<L,I)^'AN*OEtl+B*^(OEL2-DELO) 

10 CONTINUE 
XP«=^PHT( 1) 

YPM==PMT(2) 

RETURN 

901 HRITEC6,9001J 

9001 FORHAKYOHOTABLES CF POLAR MOTION COVER ONLY FROM 1958.0 TO 1970.5 
IPLEASE EXTEND.) 

STOP 

END 
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SUBROUTINE ; PRECSS 


CALL STATEMENT : PR E CSS (ET, PRE, DPRE, KK) 

SUBROUTINE PURPOSE : Calculates the precession matrix P and its time 

derivative P. P is an orthogonal matrix that transforms a vector from 
the mean equatorial system of 1950. 0 into the mean equatorial system 
of date (ET). 

INPUT PARAMETERS : 

A. Scalars. 

1. ET (extended precision). The epoch for which P and 
P are required in Julian days minus 2440000. 0, 

2. KK (integer). An index to be set to 1 when only P 
is required or to 2 if P and P are required. 

OUTPUT PARAMETERS : 

A. Matrices (extended precision). 

1. PRE (3 X 3). The precession matrix P. 

2. DPRE (3 X 3). The time derivative of the P matrix 
or P. 

PROGRAM DESCRIPTION ; Outlined in [Mueller, 1969] Pages 62-65. 

SUBROUTINES REQUIRED ; 

A. O. S. U. Project Library: 

1. ROTATE 

2. TRIM 

B. Fortran Scientific Subroutine Package: 

1. DGMPRD 

REFERENCES : 

A. Mendez, J. C. and R. J. Stern (1969). "Geographic and 
Selenodetic Coordinate Transformation Programs, " TRW N 69, 
FMT 749. 

B. Mueller, Ivan I. (1969). "Spherical and Practical Astronomy 
as Applied to Geodesy, " Frederick Ungar Publishing Co. , New 
York. 


Preceding page blank 
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PRECSS 


C 

c 

c 

c 

c 

c 


c 


SUBROUTINE PRECSS {ET,PREtDPRE,KK) 

CALCULATES PRECESSION MATRIX IN RADIANS AND ITS TIME DERIVATIVE 
IN RADIANS PER EPH6HERIS DAY ! THE PRE MATRIX PcRFCRMS TRANSFORMATION 
FROM MEAN EQUATORIAL OF 1950.0 { JED=2A332B2.423 ) TO MEAN EQUATORIAL 
OF DATE <ET) 

KX=1 IF ONLY PRECESSION MATRIX IS DESIRED 

KK=2 IF PRECESSICN MATRIX AND ITS TIME DERIVATIVE ARE DESIRED 

DIMENSION CP(3»AJ»DPt3iA) sT lA ) »PRE( 3»3 ) »OPRE( 3t31 »AH3)»DH3)»A2(3 
P,3),A3(3,3),AAI3,3),A5(3,3),At,I3,3) 

THE MATRIX CP CALCULATES THE ANCLES ZETAQ t TfcTA t Z£1 

DATA CP /3*0.U0,230A.951fe7D0,200A.2582600,250A.95162D0,.3021650U,- 
P,42668500, 1.09519500,. 018D0,-.0AieD0,.0183200/» OP ^2304. 9516700,2 

PO 04. 25826D0, 230A. 9516200,. 60A3300,-. 8537700,2. 1903900,. 05AD0,-. 125 

P4D0, .0549600,3*0.00/ ,£>1NTC /36524.21987900/ 

P,RAOS£/206264. 80624709655500/ 

TD=ET-24332B2 .42300 


TCeTD/OlNTC 
Ttl J=1.D0/RADSE 
T(2J=TC*T(1> 

TI3)=TC*TI2> 

T(4}=TC*t(3) 

CALL0GHPR0(CP»T,A1,3,4, 1) 
CALLOGHPRD (DP, 7,01,3, 4,1) 


0011=1,3 

1 Dl(I)=01(n/OiNTC 

CALLR0TATE13,-A1(3) ,A3> 
CALLR0TATE(2, A1(2),A2) 
CALL0GHPRD(A3,A2,Ah, 3,3,3» 
CALLR0TA7E(3,-A1{1) ,A5» 
CALL0GMPRD(A4,A5,PRE»3,3»3) 
I FtKK.NE .2 1G0TO7O7 
CALLLUCA(3,-01«3) ,A2) 
CALL0GMPR0(A2,PRE,A3,3,3,3) 
CALLLUCA(2,DU2),A2) 
CALLDGMPR0(A2, A5,A6,3,3,3) 

C ALLDGHPRO ( A4 , A6, AS , 3 ,3 , 3 ) 
CALLLUCA(3,-OlliJ ,A2) 
CALLDGHPR0(PRE,A2,A4,3,3,3) 
CALL0GHAD0{A3,A4,A6,3,3) 
CALLDGHAD0(A6,A5,DPRE,3,3) 
707 CONTINUE 
RETURN 
END 


2 
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SUBROUTINE : PREITR 


CALL STATEMENT : PREITR (ET, EMU, EKE, N, AM, P) 

SUBROUTINE PURPOSE : Prepares elements for STVITR subroutine for the 
Keplerian motion of a satellite about a primary body. 

INPUT PARAMETERS : 

A. ET is the initial epoch in Julian days minus 2440000. 0. 

B. EMU is the product of the gravitational constant and the 
total mass of the system. 

C. EKE is a vector of Keplerian Orbital Elements in extended 
precision and in radians or radians per day. 

1. EKE(l) is the longitude of the ascending node (Q). 

2. EKE (2) is the argument of perigee (co). 

3. EICE(3) is the inclination (i). 

4. EKE(4) is the eccentricity (e). 

5 .EKE (5) is the mean anomaly at epoch ET. 

6. EKE(6) is the mean motion (n). 

OUTPUT PARAMETERS : 

A. P is a vector of output parameters: 

1. P(l) is a precision indicator, defined as: 

-N 

P(l) = 10 

2. P(2) is the initial epoch (ET) in Julian minus 2440000. r- 

3. P(3) is the orbital major semi axis computed from: 

P<3, =[-^f 

where fi = k®(E + M), and k® is the Gaussian Gramta- 
• tional Constant. 

4. P(4 to 6) are the same as EKE (4 to 6). 

B. The AM 3x3 matrix is formed from: 

AM = R3(-*Q) Ki(-i) R3i-<^) 

AM is the transformation matrix from x'y^z' (orbital plane 
perigee) system into the KYZ ecliptic or equatorial celestial 
system. 
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PROGRAM DESCRIPTION : For a given Keplerian orbit of a satellite (ET, 
EMU, EKE) the routine sets the AM, P output parameters which are 
then transferred by the calling program to the STVITR routine every 
time the state vector of the satellite is required. 

SUBROUTINES REQUIRED : 

O. S.U. Project Library; 

• 1,. ROTATE 
2. TRIM 


REFERENCES ; 

Mueller, Ivan I. (1964). ’Introduction to Satellite Geodesy, ” 
Frederick Ungar Publishing Co. , New York. 


PREITR 


SUBROUTINE PREITR ( ET,EMU,EKE»N,AM»P ) 

IMPLICIT REAL+8 (A-H,0-Z) 

DIMENSION AMO t3) tB ( 3 t3) tC(3,3) »EKE(6) tPi 6) 

P(1)=1«D0/10.D0=*=*N 

P(2)=ET 

P {3}=(EMU/(EKE(6)*EK£t6) ) )**.3333333333333333D0 
DO 1 1=4,6 
PtI)=EKE(I) . 

CALL ROTATE (3, -EKE (2 ), AH) 

CALL ROTATE ( 1 ,-EKE ( 3 ) ,B ) 

CALL TRIM {B,AiM,C) 

CALL ROTATE ( 3 ,-EKE ( 1 > , b ) 

CALL TRIM {B,C,AM) 

RETURN 

END 
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SUBROUTINE : PRESTV 


CALL STATEMENT : PRESTV (ET, EM, EK) 


ROUTINE PURPOSE : Generation, of series of constants and a transfor- 
mation matrix for the calculation of the state’ vector of a satellite in a 
Keplerian orbit about the primary body. {The state vector is created 
by use of a companion subroutine STVKEP). 


INPUT PARAMETERS 
A. Scalar 

1 . 


2 . 


B- Vector 

1 . 


(extended precision). 

ET is the initial epoch in Julian days minus 2440000. 0. 
EM is the product of the gravitational constant and the 
total mass of the two bodies. 

(extended precision). 

EK(6) is a 6-element vector containing the Keplerian 
orbital elements in the following order and in radians 
or radians per day; 

a. Longitude of ascending node (Q). 

b. Argument of Perigee (co). 

c. Inclination (i). 

d. Eccentricity (e). 

e. Mean anomaly (M) at epoch ET. 

f. Mean motion (n). 


OUTPUT PARAMETERS 

A. Scalars 

1 . 

2 . 

3. 

4. 

5. 

6 . 

7. 

B. Vectors 

1 . 


: (All in common storage sector /EPLER/). 

(extended precision). 

E is the eccentricity. 

E2 is the eccentricity squared. 

RZ is 1 + e®''®. 

ELO is the mean anomaly (M) at the standard epoch set 
equal to Input Parameter EK(5). 

ETO is the mean motion (set equal to Input Parameter E'. 
EN is the mean motion set equal to Input Parameters EK(- 
AO is the orbit's semi major axis. 

(extended precision), 

R(l) through R(7) are coefficients of a series expansion 
for the radius vector. 
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c. 


2. F(l) through F(7) are coefficients of a series 

expansion for the true anomaly (see f Brouwer and 
Clemence, 1961]). 

Matrix (extended precision). 

1. Tl (3 X 3) is an orthogonal transformation matrix to rotate 
a vector from the x'y^z^ (satellite orbital plane) 
system to the xyz (fixed ecliptic) system. 


PROGRAM DESCRIPTION : 

A. The radius vector semi major axis ratio is given in the form 
[Brouwer and Clemence, 1961] Equation 73, Page 76: 


X 

= Ri + R3C0S‘t+ RsCOsB-t* * * 

where Ri are functions of the eccentricity and -t is the mean 
anomaly, for example: 

, 1 
Ri = l+Y 


Ra “= -e + 


3 3 _ 5 5 7 

8 192 9216 


B. The true anomaly (f) is expressed by fBrouwer and Clemence, 
1961] Equation 75, Page 77. 

f = + Fi sin + Fs sin 2-t. • • • 


where F^ are functions of the eccentricity, for example 


Fi 


2 e 


1 3 



+ 


107 7 

4608 ® . 


C. In both the R and F coefficients one must note that the eccen-' 
tricity terms are truncated beyond the e^ term. Therefore, 
for orbits with large eccentricites the effect of higher order 
terms should be closely investigated. 

D. A rotation matrix to rotate coordinates in the orbital plane system 
to the mean ecliptic system is formed in Tl as; 


Tl (3, 3) = R3(-n) Ri(-i) R3 (-w). 

SUBROUTINES REQUIRED : 

O. S.U. Project Library 

1. ROTATE 

2. TRIM 
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REFERENCES : 


Brouwer, D. and G. Clemence (1961), "Methods of 
Celestial Mechanics, " Academic Press, New York, 
pp. 76-77. 
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PHESTV 


SUBROUTINE RRESTV {&Tt£M,n<> 

PREPARING ELEMENTS FOR STVitE!* SUBROUTINE IN COMMON /EPLER/ 

£T EPOCH IN CJ0~2440000») 

Bh PRODUCT OF GRAVITATIONAL CONSTANT AND TOTAL MASS OF THE THO BODIES 
EK KEPLERIAN ORBIT ELEMENTS ir THE FOLLOWING ORDER s 

LONGITUDE OF NODE , ARGUMENT OF PERIGEE , INCLINATION 
ECCENTRICITY , MEAN ANOMALY » MEAN MOTION 
IKPL1C1TREALA8(A-H,0-I} 

DIMENSION Tl(3,3J,T2(J.3),T3I3,3J,RI7).F(7},EKt6l 
COMMON /EPLER/ E»E2tRZ»R.F,fcL0,eTU,€N,A0,Tl 

DATA ONE,THO,TRItFOR,FIFtSEV/l.OO,2.0O,3.DO,4.OO,5.DO,7.OO/ 

C SERIES FOR RADIUS VECTOR 
E=EK(4» 

E2=£*E 

E3=E2*E 

e4=E2*E2 

E5=E3*E2 

E6=E3*E3 

E7=E4*£3 

R2=ONE+E2/THO 

R(1 )sSEV*E7/9216.-FlF*E5/192.D0*TRI»E3/a.00-E 
R <2 } =~E2/TW0+E4/TRI-E6/i6 . 

RC3J=-TRI*E3/8.00-*45.00*E‘>/ 128.00-567 .AE7/5120. 

Rt4J=-E4/TRI+TH0*E6/FlF 

R{5)»— 125.D0*E5/3a4.OQ+4375.*E7/9216. 

RC6l=-27.*Eb/80. 

R{77=-16607.»£7/46080. 

C SERIES FOR TRUE ANOMALY 

F f 1 }=TH0=»E-E3/F0R+F IF^^ES/VB .00+ 107. *E 7/4608 . 
F{2J*FIF»E2/FOR-11.DO’>E4/24.00+17.*E6/192. 
FC3Jal3.DO*E3/12.DO-43.DO=!'E5/6‘».DO+95.*E7/512- 
F (41 =103 .D0+E4/96 .DO-45 1 . *Eft/4E 0 . 
Ff5)=l097.00*€5/960.00-5957.*E7/46Qa. 

F(61=l223.*E6/960. 

F(71=47273.#E7/32256. 

CALL ROTATE (3t“EKC2) tTl) 

CALL ROTATE (1,-£K{3) ,T2 J 
CALL TRIM (T2tTl,T3) 

CALL ROTATE- (3, -EK{1),T2J 
CALL TRIM (T2,T3*T1J 
EN=£K(6> 

ETO=ET 

EL0=EK{5) 

A D= ( EM/ ( EN»EN n *♦ .333333333333333300 
RETUiy< 

END 
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SUBROUTINE ; REDPI 

CALL STATEMENT : REDPI (YTH) 

SUBROUTINE PURPOSE : To reduce a given angle 0 to the interval 
0 < 0 < 27T 


INPUT PARAMETERS ; 

A. Scalar (extended precision), 

1. YTH is the input angle in radians. 

OUTPUT PARAMETERS : 

A. Scalar (extended precision). 

1, YTH is returned from the subroutine in the desired 
range (in radians). 

PROGRAM DESCRIPTION : None 


SUBROUTINES REQUIRED : None 


REFERENCES ; None 
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RE DPI 


SUBROUTINE REDPI (YTH) 

C REDUCTION OF ANGLE YTH TO THE INTERVAL ZERO AND TWO PI 

REALMS YTH*PI2 
P 12=6.28318530717958600 
YTH=YTH-PI2=!'DFL0ATf 1DINT(YTH/PI2) ) 

10 IF(Y7H .LE.O.DO)GOT013 

11 IFIYTH .LT.Pi2)G0T012 

YTH =YTH -PI2 

GOTOll 

13 YTH =YTH +PI2 

GQTOlO 

12 CONTINUt 
RETURN 
END 
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SUBROUTINE : ROTATE 


CALL STATEMENT : ROTATE (N, ANG, RNA) 

SUBROUTINE PURPOSE ; Forms a 3X3 rotation matrix, RNA, representing 
a rotation of angle ANG (in radians) about an axis xyz, designated 
1, 2 and 3 respectively. 

INPUT PARAMETERS : 

A. Axis about which rotation is to occur N. 

B. Rotation angle in radians ANG which may be positive or 
negative in sign. 

OUTPUT PARAMETERS : 

A. Rotation matrix RNA (3X3). 

PROGRAM DESCRIPTION ; 

A, Rj ( 0 ) is formed where i = N, 6 = ANG and Rj ( 6 ) = RNA. 

B. The integer variables Nl and N2 are set to the following value 
dependent on N (by use of integer arithmetic) 

N Nl N2 

12 3 

2 3 1 

3 12 

Then after RNA is zeroed : 

RNA (N, N) = 1 

RNA(N1, Nl)= cos (ANG) 

’rNA(N 2, N2) = cos (ANG) 

RNA (Nl, N2) = sin (ANG) 

RNA (N2, Nl) = -sin (ANG) 

SUBPROGRAMS REQUIRED : None 

REFERENCES : 

Mueller, I. I. (1969). "Spherical and Practical Astroi»my 
as Applied to Geodesy" Frederick Ungar Publishing Co. , 

New York, p. 43. 
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ROTATE 


SUBRCUTINER0TATt(N9ANGTRNA) 

C ROTATION MATRIX «RNA« ABOUT AXIS ”N'* 3Y ANGLE ••ANG" 

IHPLICITREAL^8IA-H,0-2) 

DIMENSI0NRNA(3,3) 

0D304Il=lf 3 
DCi303I2=lf3 

303 RNAt IltI2)=0.D0 

304 CONTINUE 
SANG=OSINUNG) 

CANGs=DCOS{ANG) 

RNA<N»N)=1*D0 

NlsN+l-( <N+l)/4)=»‘3 

N2=N+2*“({N+2)/4)=f3 

RNACNltNl)=CANG 

RNA(N2,N2)=CANG 

RNA<N1,N2>=SANG 

RNA«M2,Nl)=-SANG 

RETURN* 

END 
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SUBROUTINE : SEFODI 


CALL STATEMENT : SEFODI (M, N, K, L, Kl, K2, FF, FD) 

SUBROUTINE PURPOSE : Generation of second and fourth modified differences 
for Everett interpolation of L sets of tabulated quantities. 


INPUT PARAMETERS : 

A. Integers (single precision). 

1. M is the number of a tape unit which contains L sets 

of K tabulated ephemeris quantities written in binary mo 

2. N is the number of an output storage unit (disk) on whic • 
the function table with second and fourth differences are 
to be written in binary mode. 

3. K is the number of variables in a set. For example, ir 
the tabulated quantites are components of a state vector, 
K should be set to 6. 


4. Kl and K2 must be set to 9 and 3 respectively. 

B. Matrices (extended precision). 

1. FF (9 X K) contains 9 rows of K tabulated quantities 
input from the unit designed by M, The second and 
fourth modified differences are computed for the middle 
row (i. e. , row 5). 


OUTPUT PARAMETERS : 


A. 


Matrices (extended 
1. FD (3 X K) 

a. 


b. 

c. 


precision). 

contains the following information in rows: 
The K tabulated quantites from unit M 
(identical to row 5 of FF). 

The K modified second differences. 

The K modified fourth differences. 


PROGRAM DESCRIPTION : The second and fourth modified differences are 

computed as outlined in f O' Handley, 1969]. The method is intended to 
facilitate the use of Everett's fifth order interpolation formula (see 
EVERAL). 

SUBROUTINES REQUIRED : 

Fortran Scientific Subroutine Package: 

1. DGMPRD 
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REFERENCES ; 

O'Handley, Douglas A. et al. , (1969). "JPL Development 
Ephemeris Number 69{' JPL Technical Report 32-1465, 


SEFODI 


I, 


SUBROUTINE SEFOOI (M,N,K,LtKl»R2»FF,F0l 

GENERATING S cC -JU AND FOURTH MODIFIED DIFFERENCES htwarv 

H — IWIT CONTAINING TABULATED QUANTITIES IN L ROWj Or K EACH IN BINARY 
N - UNIT TO WRITE IN BINARY MODE THE TABLE. OF DIFFERENCES 
IMPLICIT REALMS {A-H»0-2) 

DIMENSION FF(K1,K).FD(K2,K) »CC(3,9) . 

0ATACCtl,lltCCU.2)tCC(l.3),CCU,A»tCC(l,6),CC(l»7),CC(l,b>iCCtl,9 

P),CCU,5).CC(2,1) tCCI2,9),CC(2,2l (2, 8 1 ,CCC2,3) ,CC 1 2 ,7 > ,CC 1 2,4 W 

PCC{2,6),CC{2,51,CC(3,1),CC{3,9},CC13,21,CC(3,8J,CC{3,3),CC(3,/).CC 
P(3,41,CC(3,6),CC{3,5}/e*0,DC,l.DO,2*.A2990-2,2*-.047i>12UO,2^.199L9 
P200, 2*. 56245600,-1.4366700, 2* .Ct>b489DO, 2^-. 82616ieO,2»4.B87-Ot>DO, 2 

P*-l2.0O94l9DO,lS.35961D0/ 


REWIND H 
REWIND N 
00 1 1=1,9 

1 REAOIHI {FFa,K0),K0=l,R) 
NK=L-9 


00 2 1=1, NK 

CALL0GMPRD(CC,FF,F0,3,9,>O 
DO 3 KO=l»K 
DO 3 J=l,8 

3 FF(J,R0)=FF{J+1,K.0» 

READ<H,END=9nFFl9,K0),K0=l,KJ 
2 HRITE(N) ( (FOl J.KOI , J=l ,3) ,K0=1,K.) 
9 REWIND N 
RETURN 
END 
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SUBROUTINE : SKEPTR 


CALL STATEMENT : SKEPTR (S, M, K) 

SUBROUTINE PURPOSE : Transformation of a satellite state vector into 
instantaneous Keplerian elements. 

INPUT PARAMETERS : 

A. Scalar (extended precision). 

1, M is the product of the universal gravitational constant 
(k^ and the combined masses (mi + mg) of the primary 
body and the satellite. 

B. Vector. 

1. S is a six element vector containing the state vector 
in linear units and linear units per day. The units in 
the state vector must be compatible with those in the 
scaler M, above. 


OUTPUT PARAMETERS : 

A. Vector (extended precision). 

1. K is a six element vector containing the elements of 
the instantaneous Keplerian orbit in radians/radians 
per day, in the following order: 

a. The longitude of node, 

b. The argument of perigee, 

c. The inclination. 

d. the eccentricity. 

e. The mean anomaly. 

f. The mean motion. 

2. The following factors are used in the programs (within 
the extent of IBM 360 extended precision): 

a. If the inclination is zero, then the longitude 
of node is set to zero. 

b. If the eccentricity is zero, then the argument of 
perigee is set to zero. 

PROGRAM DESCRIPTION : See references. 

SUBROUTINES REQUIRED : None. 
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REFERENCES ; 

A. Escobal, P. R. (1965). '‘Methods of Orbit Determination, “ 
John Wiley and Sons, New York. 

B, Mueller, I. I. (1964). "Introduction to Satellite Geodesy, " 
Frederick Ungar Publishing Co. , New York. 
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SKEPTH 


SUBRCUTINe SKEPTR (S,H,K) 

C TRAN^ORHATICN FROM STATE VECTOR INTO INSTANTANEOUS KEPLERIAN ELEMENTS 

C S STATE VECTOR 

C M PRODUCT OF GRAVITATIONAL CONSTANT ANO TOTAL HASS OF THE TWO BODIES 
C K KEPLERIAN ORBIT ELEMENTS IN THE FOLLOWING ORDER : 

C • LONGITUOE OF NODE , ARGUMENT OF PERIGEE » INCLINATION 

C ECCENTRICITY , MEAN ANOMALY , MEAN MOTION 

C IF INCLINATION IS ZERO - LONGITUDE OF NODE IS SET’TO ZERO 

C IF ECCENTRICITY IS ZERO - ARGUMENT OF PERIGEE IS SET TO ZERO 

C MEAN anomaly IS MEASURED FROM PERIGEE OR FROM NODE OR FROM AXIS (11 

C ARGUMENT OF PERIGEE IS MEASURED FROM NODE OR FROM AXIS U 1 

IHPLlCITREAL* 8 CA-h, 0 - 2 } 

REALMS M,K 

DIMENSION SI 6 »,Kt 6 J ,U( 3 ) »V< 3 ) 

R=O.DO 

V£= 0.00 

RO»O.DO 

DO I I=lt 3 

R=R*-S(I )**2 

VE=VE+S(l+ 31+»2 

1 R 0 ®RD+S(IJ*StI+ 3 ) 

R=DSQRTtR) 

RT=RO/R 

A=1.0O/{2.DO/R-VE/K) 

K(6)=DSQRT{M/{A-4<3n 

T=RO/OSQRTtA*Ml 

Q=I.CO-ft/A 

K (A ) =OSQRT ( T»* 2 +Q** 2 ) 

IF(K( 4 n 8 , 9,8 

8 K{ 51 =DATAN 2 «T/K{A),!)/K{ 4 ))-T 

9 P=A*(1 .D0-K{A)**2> 

DO 2 1 = 1,3 

um=sm/R 

2 VTI) = (S(l+3J*R-sni*RT)/DSQRT(P*H) 
Cl*0SQRTt{Um+V{2)l**2+{U<2>-va))**2)-1.0D 
SL=(Uf 21-VtI))/(l-UO+CIl 

CL=(uni+v(2n/(i.oo+ci) 

EL= 0 ATAN 2 CSL,CLI 
SI= 0 SQRT(U( 3 }*^ 2 +V( 3 } 4 * 2 J 
IFtSl) 3 , A , 3 

3 KC 3 )=OATAN 2 ISI,CIJ 
SU=U( 3 }/SI 
CU*Vt 3 )/SI 
UU» 0 ATAN 2 (SU,CUI 
GOTO 10 

A UU=EL 

10 IFCKtAJ) 5 , 6,5 
5 SV=RT* 0 SQRT(P/H)/MAJ 
CV=(P/R- 1 . 00 )/K(A) 

TAsOATAN 2 fSV,CV) 
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SKEPTR (Cont> 


GOTO 7 

6 TA=UU 
K(i>>=UU 

7 K(U*EL*UU 
K { 2 ) =UU-TA 

CALL REDPI (K(D) 
CALL REDPI (K(2)) 
CALL REDPI IK(Sn 
RETURN 
END 
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SUBROUTINE : SPEQU 


CALL STATEMENT : SPEQU (A, B, M, N, K, L) 

SUBROUTINE PURPOSE : Copies a two dimensional matrix (B) into layer 
L of a three dimensional matrix A. 


INPUT PARAMETERS : 

A. Scalar Integers (single precision). 

1. M, N, K. Dimensions of the three dimensional output 
matrix A (Rows, Columns, Layers). 

2. L. The number of the layer of A (consisting of K layer 5 
into which the input matrix B is to be copied. 

B. Matrices (extended precision). 

1. -B. The M X N matrix to be copied. 


OUTPUT PARAMETERS : 
A. Matrices 
1. A. 


(extended precision) . 

The output MXNXK three dimensional matrix. 


PROGRAM DESCRIPTION : None 
SUBROUTINES REQUIRED : None 
REFERENCES : None 
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SPEQU 


SUBRCUTINESPEQUU,8,M,N,K,L) 

C MAKING LAYER L Of A EQUAL TO B 

1MPHCITREAL*8 CA“H,0-*Z> 
DIMENS.I0NA(M,N,K)7B(M,N) 

OD2J=lfN 
2 A(ly 
I CONTINUE 
RETURN 
END . 
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SUBROUTINE : STRANS 


CALL STATEMENT ; STRANS (JB, GAST, SMAT) 

SUBROUTINE PURPOSE : The subroutine computes the rotation matrix S to 

rotate coordinates from the true celestial system to the average terrestrial 
system. 

INPUT PARAMETERS : 

A. Scalars (extended precision). 

1. JD is the epoch in Julian days. 

2. GAST is the Greenwich apparent siderial time in radian 

OUTPUT PARAMETERS : 

A, Matrix (extended precision). 

1. SMAT is a 3 X3 array containing the rotation matrix 
•from the true celestial system to the average terrestri.^: 
system. 

PROGRAM DESCRIPTION ; The SMAT matrix (in Mueller's notation the S matrixl 
is computed from the polar motion components (x, y) and the Greenwich 
apparent siderial time (0) from 

S = Ri(-x) Ri(-y) R3(0) — 

SUBROUTINES REQUIRED : 

O. S.U. Project Library: 

1. POLE 

2. ROTATE 

3. TRIM 


REFERENCES ; 

Mueller, I. I. (1969), "Spherical and Practical Astronomy 
as Applied to Geodes3% " Frederick Ungar Publishing Co. , 
New York. 
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STRANS 


SUBROUTINE STRANS (JD, CAST, SMAT) 
lNPLlCI7REAL=^t (A-H,0-Z) 

DOUBLE PRECISION MJD,JD 

DIMENSION R3TETA(3,3) »RlMY(3t3) ,R2MX(3*3) tTU3»3) ,SMAT(3»3) 
DATARADSE/206264 ,606247096400/ 

COMPUTES THE SMaTRIX TO ROTATE FRM TRUE CEL, TO AVERAGE TERR. SYST 
INPUT PARAM. — JDtGAST. 

EXTERNAL ROUTINES — POLE (FOR POLAR MOTION), AND TRIM 

MJO=JD”-2400000.500 

CALLPCLE(MJD,XPP,YPP) 

XPP=XPP/RAOSE 
YPP=YPP/RADSE 
C ALLROTATE ( 3,GAST,R3TETA ) 

CALLR0TATE(1,“YPP,R1MY) 

CALLROTATE (2,-XPP,R2MX) 

CALLTRIM (R2MX,RlHY,Ti) 

CALLTRIM(T1,R3TETA,SHAT) 

RETURN 

'60 FCRHAT{//,-20X,2F20*9) 

END 


- 140 - 



SUBROUTINE : STVITR 


CALL STATEMENT : STVITR (ET, STV, AM, P) 

SUBROUTINE PURPOSE : To calculate a state vector of a satellite in a 
Keplerian orbit for a given epoch and in components of a fixed 
Cartesian coordinate system. 

INPUT PARAMETERS : {Extended precision. ) 

A. ET is the epoch in Julian days (t) for which the state vector 
is required minus 24400000. 

B. AM and P are the 3x3 matrix and the six element vector 
described in PREITR. 

OUTPUT PARAMETERS : (Extended precision. ) 

A. STV is a 6-element state vector of epoch in the XYZ system 
of a satellite moving in a Keplerian orbit about a primary 
bodj\ 

PROGRAM.DESCRIPTION : . 

A. The mean anomaly (-t) for the epoch of observation is 
obtained from 

I = la + n(t - to ) 
where 0 ^ -t. < 2tr 

B. The eccentric anomaly E is then found from the solution of 
Kepler's equation by: 

1. Estimating the initial value as 

Eo = -t 

2. Then using a Newton iterative procedure to find the 
E p of epoch as follows : 

a. Compute AEj from: 

AEj = Ej -e sin El - I 
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b. Compare AEj with the precision estimator 
P ( 1 ) which is denoted P 

i AEiI < P 

c. Then: 

If I AEj 1 < P, terminate iteration and 
continue program. 

If I ilEj I s p, continue iteration by setting 

If convergence is not achieved after 50 iter- 
ations, print error message and stop. 

C. Compute the auxiliaries (notation from Esc obal). 



= 

sin Ep 

SE 

c. 


-c Computer 

cos Ep / 

Variables 

CE 

aE 


na 

CON 


1 - eC, 

Compute the 

state vector in the x', y', z' 

system: 

x' 

= 

a(Ce - e) 

1 


y' 


a(l-e^)^ S, 

XOR Vector 


- 

0 


x' 

= 

- a E • S* 


r 

- 

o A 

aE(l-e®)^ C, 

VOR Vector 

Of 

Z 


0 



E. Rotate the above state vector into the XYZ system by (see 
PREITR). 


STV 


AM 

0 


x' 


AM 

0 


XOR 



• 

y' 






0 

AM 


z' 


0 

AM 


VOR_ 


_ 


k' 






y 






7/ 
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SUBROUTINES REQUIRED : 

O. S.U. Project Library: 
1, ONEM 


REFERENCES : 

A. Escobal, P. R. (1965). "Methods of Orbit Determination, ' 
John Wiley and Sons, New York. 

B. Mueller, Ivan I. (1964). "Introduction to Satellite Geodesy 
Frederick Ungar Publishing Co. , New York. 
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STVITR 


SU3RCUTINE STVITR ( ET ,STV ,AM,P) 

IMPLICIT REAL^i^eCA-HtO-Z) 

DIMENSION AM (3,3) ,STV ( 6 ) ,XOR( 5 ) , VOR (3 ) , VI { 3 ) , V2( 3 J ,P ( 6 ) 
DATA ONE,XOR(3) ,VOR(3i/l.DO,2=4=G.DC/ 

EM=UMOD( (P C5)+P(6)*(ET-P (2) ) ) ,6 .2831S33071795&6D0 ) 

E=EM 


00 I Isl,50 

D£= ( E*-P (4)^0SIN (E )~EM) / { QNE-4^ (4>t'0C0S {E ) ) 
1F(DABS(DE).LT.P(1)} GOTO 2 

1 E=E-DE 
WRITE(6,70)DE,E 

2 SE=OSIN(E) 

CE=OCOS(E) 

CON= P{6)*P(3) /(0NE-“P(4)%CE> 

XQR( 1)=PC3)=J'(CE-P{4>) 

X0R(2J=P(3)*DSGRT{0NE“P(4)=J:P{4) )*SE 
VOR C i)=-CGiN*SE 

V0R(2)=DSgRT{0NE-P(4)=}cp(4) )«C£=i'CON 
CALL ONEM (AM,X0R,V1) 

CALL ONEM iAH,V0R,V2) 

00 3 1=1,5 

sTvm=vi{ij 

3 STV( I+3)=V2(I) 

RETURN 

70 ECR.MATdOX, »ND CONVERGENCE AFTER i?0 ITERATIONS* ,2024.16) 
END 
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SUBROUTINE: STVKEP 


CALL STATEMENT : STVKEP (ET, STV) 

SUBROUTINE PURPOSE : To calculate a state vector in a Keplerian orbit. 

(This subroutine must be used in conjunction with subroutine PRESTV 
which generates a series of coefficients in common storage used in 
the state vector computation.) 

INPUT PARAMETERS ; (Extended precision. ) 

A. Scalar. 

1. ET is the epoch for which the state vector is required 
in Julian days minus 2440000. 0. 

OUTPUT PARAMETERS : 

A. Vector (extended precision). 

1. STV ( 6 ) is the output state vector. 

COMMON AREA PARAMETERS : 

A. Area /EPLER/ E, E2, RZ, R, F, ELO, ETO, EN, AO, T3. 
The common area parameters have been described in PRESTV. 

PROGRAM DESCRIPTION : 

- A. The mean anomaly (-t) at epoch t (variable_ET) is computed 
from the mean anomaly at epoch to, (which is variable ETO 
in common storage), and n which is variable EN in common 
storage: 

E L “ ” 'Co n(t to) 

where 0 < j-tl < 2it . 

B. The S vector is then computed by: 

Si. = sin-C 

Ss = 2 sin®-t = sin 2 -t 
^ = (3 - 4sin®'C)sin-t = 

84 = (2-4 sin® V)2 sin-t 


= sin3-t 
cos-t = sin4't 


&7 = sin7-t . 
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C. The C vector is then computed from: 

Cl = cos-t 

Ca = (2 -• 4 = 1 - 2 sin^-t = cos 2 -t 


C 7 = cos7-t . 

D. The radius vector (-|-) is computed from TBrouwer and 
Clemence, 1961] Equation 73, Page 76: 

-T ' i«'°‘ • 

1=1 

E. The true anomaly (f) is computed from fBrouwer and Clemence, 
19611 Equation 75, Page 77: 

*7 

f = ^ El Si . 

F. Note in the above expressions that the H and F coefficients 
exclude terms greater than , thus the effect of the exclu- 
sion should be investigated before using this subroutine. 

G. The radius (r) is then computed from; 

r 

r • a . 

a 

H. The Cartesian coordinates in the orbital plane are then computed 
from: 


X 


rcosf 

y 

= 

r sinf 


= 

. 0 


I, The time derivatives are then formed from: 


» 


r V 1 

X 


-na -J- 

rv 1 -e 

• 


X + e a a/ 1 - e^ 

y 


n a. 

r 

z 


0 
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J. The coordinates and time derivations are then rotated into 
the XYZ system through the matrix T3 (from PRESTV) 
which is effectively: 


and 


T3 = Ra Ri (“i) Ra (-W) 
’T3 O' 


S = 


“T3 01 Tx“ 

.0 T3J Lv_ 


where S is the output vector STV. 


SUBROUTINES REQUIRED : 

0,S.U. Project Library: 

1. REDPI 

2. ONEM 


REFERENCES : 

Brouwer, D. and G. Clemence (1961). "Methods of 
Celestial Mechanics, " Academic Press, New York, 
pp, 76-77. 
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STVKKP 


SUBROUTIME STVKEP (ET,STV) 

CALCULATION OF STATE VECTOR IN A KEPLERIAN ORBIT 
1HPUCITREAL*8(A-H,0--Z) 

DIMENSION R(7)tF{7) tS(7},C{7),T3t3t3)iX0R(3l,V0Rl3)tP0SI3J,T2(3l,S 
PTV{6> 

COMMON /EPLER/ £,£2tRZ»R»F» ELO.ETO tENt A0,T3 

DATA ONE»THOtTRItFOR,FIF,SEV/1.00,2.OO,3.0O»A.O0t5.00,7.D0/ 
MULTIPLE SINES AND COSINES OF MEAN ANOMALY 
ELsELO+EN* ( ET-ETO ) 

CALL REDPl lELJ 
SL=DS1N{EL) 

CL=DCOS(EL) 

SLQ=FOR»SL*SL 

CLQ=FOR-SLQ 

SCL=TWO*SL=»CL 

S{II=SL 

St2)=SCL 

S{3)=SL=i=(TRl-SLa) 

S CA ) =SCL* ( TWO-SLQ ) 

S(5)=SL»{FIF-SLQ=*(F1F-SLQ)) 

S(6)=SCL*(TRI-CLQ*(F0R-CLQ)) 

S(7)=SL=i‘{SEV-SLC*(lA.D0-SL0»(SEV-SLQ) )1 

cm=cL 

C 12 )=. 500’*= (TWO-SLQ ) 

C(3)=CL*(CLQ-TR11 
C ( A ) =ONE-CL0*, 5D0=» ( FOR-CLO ) 

C<5)=CL=»(FIF-CLQ*{F1F-CLQ) ) 

C (6 ) =CLC*. SOO^I SEV+TWO-CLQ“» t FIF+DNE-CLQ J) -ONE 
Cl7J=CL*(CL0^{lA.C0-CLQ*tS£V-CLQ) )-SEV) 

CALCULATION OF RADIUS VECTOR AND TRUE ANOMALY 

RV=RZ 

TA=EL 

00 T 1=1.7 

RV=RV+R{1)*C(T1 

1 TA=TA*F(1 J*S(I) 

RV=RV=*AO 

X0R{1»=DC0S(TAJ*RV 

XCR(2)=DSIN(TA)*RV 

X0R{3)=0.00 

VORU)=-EN*AO#XOR{2)/(RV=»OSQR7(ONE-E2n 
VORt2)=EN«AO*(XDR(l ) +E«AO )’*=DSQRT{UNE-E2)/RV 
VOR(3)=0.00 . 
call DNEM (T3.X0R.P0S) 

CALL ONEH {T3.V0R.T2) 

DO 2 1=1.3 
J=I + 3 

STV(n=PQS(IJ 

2 STV(J)=T2(I) 

RETURN 

END 



SUBROUTINE : TRIK 


CALL STATEMENT : TRIK {Y, YP, TETA, FETA) 

SUBROUTINE PURPOSE : The subroutine is used with either FUNPL5 or 
FUNPL6 and provides the second time derivatives of the physical 
libration parameters and the ©» *1? partial derivative matrices 
(see rPapo, 1971] Pages 99-106). 

INPUT PARAMETERS : 

A. Vector (extended precision). 

1. Y is a 6-element vector containing the Eulerian angle ) 
of the Moon and their time derivatives as input through 
FUNPL5 or FUNPL6. 

OUTPUT PARAMETERS : 

A. Vector (extended precision). 

1. YP is a 3-element vector containing the second time 
derivatives of the physical libration angles 
(X, 6. p). 

B. Matrices (extended precision). 

1. TETA is a 6x6 matrix which is the© partial derivative 
array evaluated at ET fPapo, 19711. 

2. FETA is a 6x3 matrix which is-the^» partial derivative 
array evaluated at ET TPapo, 1971], 

COMMON AREA PARAMETERS : Common areas are used extensively to transfer 
data in and out of this subroutine, the calling subroutine (FUNPL5 or 
FUNPL6) and the main program. 

A. Area /MAIF/ALF, BET, GAM, TEQ, CCE, CCS, P, G. 

1. ALF, BET, GAM, TEQ are described in subroutine 
documentation for FUNPL5. 

2. . CCE is an extended precision scalar which contains the 

gravitational constant of the earth. . 

3. CCS is an extended precision scalar which contains the 
gravitational constant of the sun. 

4. P is an extended precision 3-element vector containing tt 
rotational velocities e^, By, respectively of the 
ecliptic coordinate system defined in •'Papo, 1971] 

Page 83. 
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5. G is an extended precision 3 x 3 matrix containing 
partial derivatives of r, cr , p , with respect to 
Css, , Cao respectively as defined in rpapo, 

1971] tquation 3.32.12, Page 101 . 

B. Area/EXE/ET, B, SUN. 

1. The above parameters correspond to the parameters 
TO, B, SUN described in the FUNPL5 subroutine 
documentation. 

C. Area /ALL/SEV, ESH, NUT, KLU. 

1. The above parameters correspond to the parameters 

SEV, A, NUT, KLU described in the FUNPL5 subroutine 
documentation. 

PROGRAM MATHEMATICS : The expressions used to evaluate the second partial 

derivatives are given in [papo, 1971] Pages 99-101. Comments are included 
in the statement listings to refer tfte next segment of statements to specific 
equations given in f Papo, 1971]. 

SUBROUTINES REQUIRED ; 

O, S. U. Project Library: 

1. ROTATE 

2. TRIM 

3. ONEM 

4. LUCA 

REFERENCES : 

Papo, Haim B. (1971), "Optimal Selenodetic Control, " The Ohio 
State University, Department of Geodetic Science. Report No. 156. 
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SUBROUTINE TRIK (Y,YP ,TETA,FETA ) 

SUBROUTINE TRIK FOR CALCULATING EQUATIONS OF NOTION OF THE PHYSICAL 
LIBRATION ANGLES ANO THE PARTIALS FOR GENERATING STATE TRANSITION 
AND PARAMETER SENSITIVITY MATRICES 
PROGRAMMING SECTION 3.8 IN THE DISSERTATION (PAPO) 

IMPLICIT REAL*8 (A-H,Q-Z) 

DIMENSION Y(6),YPt3) ,SEV <6 J jP < 3 1 1 B (3 J t SUN (3 J , EH IIN 1 3, 3) »0 (3) , 

PROP (3) ,OHE(3)»PRE{3),Tl{3,3) ,T213,3) ,T3(3»3) ,T4{3;3 ) ,TEM1{3I , 

PTEH2I3) ,TEH3(3)tT 5(3,3) tT6(3t3) tT7(3t3),T8(3,3)tSUN0(3) iOOO{3.3) 

P,000(3t3} tTETAISrO) >A(3) <E(3 t3> tAVI C 3t3) tRAC{3t3) 
PtGt3t3J,B0(3t3>,FETA(6,3) 

COHHDN/HAIF/ALF, BET, GAM, TEQ,CCE,CCSfPtG/EXE/ET»B, SUN/ALL/ SEV,ESH 
P,NUT,KLU . 

DATA EKiIN(l,3),£HlIN(2,3),EHllN«3,3) tPRE (3), 'ONE/1. 00, 3*0.00, 1.00/ 

Y - EULERIAN ANGLES OF MOON 

SEV - ti) LONG (2) OMEGA (3> LONG DOT C4) OMEGA DOT (5) LONG DDOT 16J OMOD 
G - MATRIX OF PARTIALS ALF,BET,GAM / C22,BET,C20 
TEQ - MEAN INCLINATION OF MOON tCUATOR 
P - ECLIPTIC ROTATION RATE VECTOR 

B,SUN - EARTH ANO SUN POSITION IN ECLIPTIC SYSTEM 
CALLR0TATE{3,YU) ,T1) 

CALLROTATE 11 ,-Y C3 ) ,T2 ) 

CALLRDTATE{3,Y(2) ,T3) 

CALL TRIM (T1,T2,TAJ . 

CALL TRIM (TA,T3,T5) 

T.S - TRANSFORMATION MATRIX H 
STEIN=0Nc/T2<3,2) 

STEIN = 1 / SIN <TETA) 

EH1IN(1,1)=T1{I,2)*T2(3,3)*STE1N 
EHlIN(l,2)=Tl(l,l)*T2(3,31*STEi;j 
EH11N(2,1)=-T1{I,2J*STEIN 
EH1IN(2,2)=-T1U,1J*STEIN 
EH1IN(3,1)=-TI(1,1) 

EH1IN{3,2J=THI,2) 

CREATING W(-l) 

0 : NEGATIVE COORDINATES OF EARTH IN TRUE SELENOGRAPHIC 
CALL ONEH (T5,B,0) 

C SUND : COORDINATES OF SUN IN TRUE SELENOGRAPHIC 

CALL ONEM (T5 , SUN , SUND) 

2E=0.D0 
ZS=O.DO 
DO 1 1=1,3 - 

00 6 J=l,6 

6 TETA(I,J)=0.00 
00 7 J=l,3 

7 FETA(I,J)=0.D0 
ZS=ZS+SUNO ( r ) *SUNO ( I J 

1 ZE=ZE+DCn*D{I J 
TETA{I,A)=1.00 
TETAt2,5»=1.00 
TETA(3,6)=1.D0 
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c 


c 

c 


c 

c 


c 


c 


c 


c 


TBIK fCotit^ 

FORHATION OF TRIVIAL PARTS IN CAPITAL FI CFETAJ ANO TETA 

Ze=CCE/ZE=»»a.5D0 

ZS=CCS/ZS**2.500 

CALCULATING THE EFFECT OF THE MOTION OF THE ECLIPTIC 
R=V(2) 

DR=y(5» 

CR=DCOS(R) 

SR*DS1N(R) 

P0Pm=-Y(6)+CR*P(l)+SR»P t2) 

PDPI2)= -SR*P(1)+CR*P(2) 

P0P(3)*Y(5J+P(3) 

POP - IS PtO) 

CALL ONEH (T4«P0P,0HE) 

0«Et3)s=0H£13J+Y(4) 

OME ~ RATES OF ROTATION ABOUT SELENOCRAPHIC AXES 
PREei)=DR*(-SR*P(l)+CR*Pt2)) 

PRE{2)=0R*(-CR*H{1)-SR^P<2)) 

PRE - PSIOOT « 0R3(PSI) / OPSI * P 
CALL ONEM (TA.PRE tTEHI) 

TEMl IS T(C) 

CALL LUCA (3,Y(4),T6) 

CALL TRIM (T6tT4,T7) 

CALL LUCA (1,-Y(6»,T6) 

T6 IS TtBJ 

CALL TRIM (T4,T6,T0> 

CALL DGMADO ( T7»T8 , T6 t 3 »3 ) 

CALL ONEH (T6,P0?,TEM2) ' 

TEH2 IS TO) * P(0) 

A ll)= 2E*D(2)*D{ BJ+IS^SUNOIZJ^SUNDOl-OMErZloOHEf 3) 

A (2)=-(ZE*D( l)*0(3)+ZS»SUNDtl)»SUN0l2 )-DM=(l)*CM£{3) ) 

A (3)= ZE*D(U*u<2)+ZS=»SUND(l)=^SUND{2)-OME(l)=»‘C;ME(2) 

A IS TtA) 

TEM3C1)=ALF*A{1)-TEM2{1)-TEH1{1) 

TEH3 (2 )=8£T*A( 2)-TEH2 12 ) -TEMl (2 ) 

T EM3 {3 ) =GAH-*‘ A ( 3 ) -TEH2 (3 ) -TcHl ( 3 ) 

TEH3 IS T(D) 

IFJKLU.GT.l) GOTO 5 

IF NO GENERATION OF PARTI ALS IS NEEDED KLU IS GT. THAN 1 

Gil=CME{2) 

G21=-OME(l) 

G12=TEHm)/OR 

G22=TEMlt2)/DR 

G32=TEMl{5)/UR 

G13=T4U,3)*P0P(2)-T4tl ,2)*POPJ3) 
G23=T4(2»3)=»P0P(2)-T4(2,2)»P0PO) 
G33=T4(3,3)*P0P(2)-T4(3t2)=»P0P(3) 

G123 ARE Gll) ROW VECTORS 

DOOt l,I)=-QKE(2)»ONE+ALF) 

DDD(2tl)- OMEd )*(ONE+3ET) 

D00(3,l)s0.00 
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00D{ lt21=-Al.F*{TA<2f3)*0H6{3J+T^I3,3l*ClMe(2J )-G12-76a,3J 
00012,2}= BET*a4}lT3)=»C.M£(3)+74(3,3)<=GM6a} )-G22-T6(2,3} 
000(3,2J=~GAM*IT4 {1,3I*0H£(2)+T4<2,3)*CH£(U J-G32-T6 (3,3) 

ODDI 1,3 )= ALF*IT4(2tnvQME(3)+T4(3, 1 )=»0.ME(2) J-G13+T6 ( 1 ,1 ) 
DOOt2,3)=-BET*(T4«l,l)*OK={3)+T4(3,l)*OKeU)}-G23+T&12,l) 

ODD 13, 3 )= CAH*(T4(l,l)»0ME{2)+T4t2,l)*0HE(lJ>-G33+Tot3»l) 

DOD IS DTIDJ / 0 EUURIAN DOTED 

CALL TRIM (E«11N,DDD,AVI> 

AVI ARE PARTIALS OF EUUERIAN DOUBLE DOTS VS. EULERIAN DOTS 

C11=2E*(T5(2,1)*B U)+T5(2,2)=B I21+T5I2,3)=*B (3)) 

C21=2E»{-T5U.i>»S U)-T5(l,2)*B (2)-T5U,3)*8 (3)) 

C12=2Et|-T5a,2)'>6 m+T5(l,H*B <2)1 
C22=ZE*(-T5(2,2)*B m+T5(2,ll=B (2>) 

C32=26’>'(-TS(3,2)»6 U)+TJ>(3, U’^B (2J) 

C13=2E=*lT4(l,3J*T3{2,l)>>6{U+r4U,3)*T3{2,2J4B{2}-T4a,2J=&l3)) 
C23=2E*CT4{2,3)*T3(2,1)=B(U+T4(2,3)*T3C2,2)*B{2)-T4{2,2)*3(3)1 
C33=2E*IT413,314T3<2,1)*8{1)+T4{3,3)*T3{2,2)*BC2)-T4I3,2»*S13)J 
C123 ARE C(l) ROMS OF EARTH 

Hll=2S*{T5(2,l)=i'SUNaj+T5(2,23*SyN(2J-»T5{2,3)*SUN|3)) 
H21»2$*(-T3(l,l) + SUN(l)-T5(i,2J»SUNI2)-T5a,3)*SUNl3>) 
H12=ZS=»{-T5C1,2}*SUN(1)+T5U,1)«-SUNI2)} 
H22=ZS*(-T512t2)*5UNll)+T5(2,i)*SUNI2)> 
H32=2S»<-T5(3,2>^5UNU)+T5(3»IJ*SUN(2)) 

Hi3=2S*(T4n,3}*T3l2,l)«SUN(l)+T4U,3}*T3l2,2)*SUNI2)-T4{l,2)4SUM 

P(3}) 

H23=2S*IT4{2,3J*T3(2,l)=SUNU)+T4£2,3}«T3l2,2)*SUN{2J-T4C2,2l*SUN 

PI3)) 

H33=2S*(T4l3,3}*T3(2,l)*SUWm+T4f3,3)*T3l2,2)«SUNI2J-T4C3,2)*SUN 

PI3J) 

H123 ARE CII> ROWS OF SUN 

00011.1) '= ALF*(C2l^0t3)+H21*SUN0l3)-G2l^0HEl3J)-TEH2t2)-^’Hif2) 
000(2,i)=-flET*(Cil*0(3)+HXl*SUN0(3)-GH*0M6|-3)l+TEM2(l)-*-TEHlti) 

00013.1) = GAH*(C11*D<2)+C21*ua J+Hn*SUNDI2)+H21*SUNOU)-GH=0HEI2 
P)-G21*0MEU)) 

OOOt 1,2)= ALF=»(C22t‘0(3)+C32»DC2)+H22»SUNO(3l+H32»SUN0<2)-G22*OM6t 
P3 )-G32*OME(2) )+T4(l,2)*PREm-T4(l,l)=fPREC2)-{T6(l,ll*PREm 
P+T6I l,2)=PRE(2 JI/OR 

0G0(2,21=-3ET=(C12=D(3)+C32^0(1)+H12*S0NDI3}+H32*SUND{U-G12^0ME 
PC3)-G32*CHE(l) !»T4(2,2)*PR£UJ-r4l2,l)«PR£t2J-{T6(2,l)*PRfcU) 

P+T6 1 2, 2 )*PR6 1 2 U /DR 

DOD (3,2)= GAH={CX2=D(2)+C22*D(l)+H12*SUND(2)+H22*SUNO(l)-Gl2=OME 
P(2)-G22*0MEa) ) ♦T-‘-(3,2)*PR£tl)-T4(3,l)“»PRe(2>-CTG(3, U*PR£IU 
P+T6(3,2)*PR£(2))/tR 

00D(l,3)= ALF^(C23»ai3)+t33'»U(2l+H23=SUNOI3)+H33=»SUNl}(2}-G23»GMfc(3 
P)-G33*0HE(2 J)-T4( l,3)*PRE (2 )+T4( 1 ,2 I*PR6f 3) +T6( 1,2)*PCP (3)-T6( i,3) 
P*P0P(2) 

DOO ( 2, 3 )=-6ET* (Cl 3*D( 3) +C33»0 ( 1 ) +H13*SUN0 ( 3 ) +H33*SUN0( I ) -G 13*0Ht (3 
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no non 


TRIK fCont) 


PJ-G33*OMEm)-T4(2i3)*PRE(2)4T4(2,2)*PRet3)+76{2t2)*P0P«3)-T6t2i3J 

P>(tpOPI2) 

OODl 3f 3 ) = GAH* (C13=*0 « 2 ) +C23+G ( 1 » +H13*St’ND ( 2 > +H23#SUN0 ( 1 ) -G13*0Mfc 1 2 
PJ-G23*0MEC1J )-TA{3,3}*PRE(2HT4(3 t 2)*PRE(3)+T6(3»2J*P0P{3)-T6(3»3) 
P*P0Pf2» 

• DOD ARE PARTIALS OF T(D) VS. EULERIAN ANGLES 
CALL TRIM (EHllNt000,T7) 

RACU»l) = (TiaTlJ*T2(3»3)»TEH3{l)-Tl{l,2}*T2C3*3)*'TEH3(2))*STtlN 
P+T7(l»l) 

RACC 2 ,l) = (-Tmtl)*TEH3m + Tl(1.2)»TEH3t2))=J'STcIN*T7{2,l) 

RAC(3,1»=T1(1»2)+TEM3(1) +TX tl , 1 )*TEH3 (2 J +T7«3, 1 ) 

RACil,2)=T7{l,2) . ■ 

RAC12,2>=T7{2,2) 

RAC(3,2)=T7(3t2) 

RAC* 1,3» = (-TU1,2)*TEH3I1 J-TlIl.D^TEHBta) )*STE1N*STEIN+T7{1»3] 
RAC{2r31=(Tl(l,2)*T2{3,3)'»TEH3(l)+Tl(l,lI*T2(3,3)*TEH3(2)) 
•P*STEIN*STEIN+T7<2t3l 
RAC(3,3)=TT{3,3) 

RAC ARE PARTIALS OF EULERIAN DOUBLE DOTS VS. EULERIAN ANGLES 
T£TAC4»l)=RACU»l)+RAC«2f 1) 

TETA<4,2)=RAC(l,2)+RAC(2,2)-TETA(A,l) 

.T£TA{4,3)=RACU,3)+RAC(2»3) 

TETA(A,4)=AVI(1,1I+AV1(2,1) 

TETAC4,5)=AVI11»2)+AVI(2,2»-TETA(4,4) 

TETA*A,6)=AVI(l,3)+AVn2,3) 

TETA(5,1)=RAC(2.1) 

TETAt5»2)=RACt 2»2)-RACt2,lI 
TETA*5t3)=RAC(2t3) 

TETAI5,4>=AVH2,1) 

T£TAC5,5)=AVI{2,2}-AVI(2.1J 

TETA(5»6)=AVI(2,3) 

TETA(6,1)=RAC(3,1) 

TETA{6,2)=RAC(3f2J-RAC(3»l) 

TETA{6,3)=RAC(3,3) 

TETA*6,4)=AVH3,1) 

TETA(6,5)=AVI(3»2)-AVl(3,l) 

T£TA*6,6)=AVI{3,3) 

FILLING IN THE 4t5f6 ROWS OF CAPITAL TETA 

00 3 1=1,3 
00 2 J=l,3 

2 EtJ,I)=EMlIN(J,IJ*A<n 

3 EC1,I)=EU,1)+E{2,I } 

CALL TRIM (E,G,BOJ 
00 4 1=1,3 

DO 4 J=l,3 
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TRIK fContt 


FETAlI-»3f J)=SO(lf J) 

FILLING IN THE 4»5»6 ROMS OF FETA 
CALL ONEM {EHlINi TEM3,TEH1) 

YPU )=TEM1(1)-S£VI5>+TEMH2) 

YP123=T£Hl(21-StV{6J 

YPI3)=TEM1(3) 

YP - ENCKE EQUATIONS OF MOTION OF 

RETURN 

END 


CCAPITAL Fl> 


PHYSICAL LIBRATION ANGLES 



SUBROUTINE ; TRIM 


CALL STATEMENT : TRIM (A, B, C) 

SUBROUTINE PURPOSE : Multiplies two 3x3 matrices, A and B, to form 
a 3 x 3 product C. 

INPUT PARAMETERS : 

A. Matrix A (3, 3) 

B. Matrix B ( 3 , 3 ) 

OUTPUT PARAMETERS : 

A. Matrix C (3, 3) 

PROGRAM DESCRIPTION : 

A. 3O3 “ 3^3 3B3 

SUBROUTINES REQUIRED : None 
REFERENCES : None 


Ptecefling page Wank 
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TRIM 


SUBRC3UTINE TKIM (A,B,C) 

REAL’^8A(3,3),Cv(3,3) ,C{5,3) 

00lI=Xf3 

D01J=1,3 

1 Ca, J)=A(I,l)«B(l?J)+AIl,2)=f=B{2f J)+A(l ,3)=i-?i3tJ) 
RETURN 
END 
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DESCRIPTION OF MAIN PROGRAMS 
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MAIN PROGRAM INDEX 


PROGRAM PAGE 

A. Simulated Ephemeris Generator 163 

B. Optical Observation Generator 169 

C. Adjustment of Lunar Phj^sical Libration Angles 175 

D. Optical Observation Adjustment 181 

E. Lunar Range Simulator 187 

F. Lunar Range Adjustment 193 


Preceding page Wank 
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MAIN PROG AM : A 


PROGRAM PURPOSE : Generation of the simulated ephemeris of the earth- 
moon system described in [Papo, 19711. 

INPUT PARAMETERS : An example of card input parameters is given in 
Appendix A/2. 

A. Vectors (extended precision). 

1. PAR is a 5-element vector containing: 

a. The initial epoch in Julian days minus 
2440000. 0. 

b. The final epoch is Julian days minus 
2440000. 0. 

c. The initial step size for the integration in 
days. 

d. Elements 4 and 5 are not used, 

B. Scalars (single precision). 

1. HMIN is the minimum step size to be allowed in the 
integration routine. 

2. UMAX is -the maximum step size to be allowed in the 
integration routine. 

3. RATI is the relative accuracy required (number of 
correct significant figures) for the geocentric state 
vector ot the moon. 

4. RAT2 is the relative accuracy required for the 
numerical integration of the earth and moon Eulerian 
angles. 

OUTPUT PARAMETERS t 

A. The main output from the routine is the state vector of the 

earth-moon system, the Eulerian angles of the system and time 
derivatives of the vectors and angles. The units and sequence 
of the quantities are described in the subroutine FUNEPH descrip- 
tion (for the Y and YP vectors) and the subroutine OUTEPH 
(for the FD matrix layered into the FH array). In addition the 
FH array of the simulated ephemeris is written on a direct 
access disk device identified in the program as unit 4. 

PROGRAM DESCRIPTION : A flow chart of the program logic is given in 
Appendix A/l. 


Preceding page blank 
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SUBROUTINES REQUIRED : 

O.S.U. Project library: 

1. EUNEPH 

2. EKHARD 
3» OUTEPH 

4. PREITR 

5. DVDPFI 

6. MEANAN 

7. RE DPI 

8. SEFODI 

REFERENCES ; ^ „ 

Papo, Haim B. (1971). ’’Optimal Selenodetic Control, 

The Ohio State University, Department oE Geodetie Science , 

Report No. 156. 
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Program A 


GENERATION OF SIMULATED EPHEHERIS OF THE MOON+EARTH SYSTEM 
IMPLICITREAL*S(A-HfO-Z) 

REAL*4 £F,HM1N,HMAX,RAT1,RAT2 

DIMENSION YC16) ,YP(18),KQ(16)«0T{20,1S) ♦ENCKEC7J .EKH(6),AttH(3»3) 
P,PM(6>,eP(18>,YNN(10) ,K0U8),PARI5)fPLA{31tDNM(6),ANH(6) 

DIMENSION FF{9,18),F0(3»18) ,FH(3,18,15) 

COMMON/MA1FUN/EMCKE,AHM»PH 
EXTERNAL FUNEPH.OUTEPH 

DATAZ£RO,PI t EKM{3 ) * EKMt Al/0 ,D0»3. 1A159265358979300 ,.08980D0t .0549D 
PO/ 

REAO(S,Al) PAR 

READ (5*61) HMI»N,HHAXfDELTTRATltRAT2 
DO:cPAR(l) 

E JD=OO+2WO0CO .DO 

CET= tDO+24980 . DO) Z36525. DO 

CALL MEANAN {EJD,ANH,OMH) 

ENCKE(1>=00 

ENCKE(2)=PI+ANMt2)-ANH(51 
CALL REDPI (ENCKEC2)) 

ENCKEt3)=ANMC5J 

ENCK6<4)=.23D0 

ENCKEJ 5 » = { OHOD ( OOf 1 .00)*2.DC+ C23925.83600+8640184.542DO*CE7+.09290 
P0*CET»CETJ/432O0.00)*PI 
CALL REDPI (ENCKE(5)1 
ENCKE<6)=. 4091600 
ENCKE(7)=6.300388098DO 
EKKI1)=ANH{5) 

CALL REDPI (EKM(D) 

EKM(2)=ANM(4)-ANH{i>) 

CALL REDPI (EKH(2)) 

EKHC5)=ANHt2)-ANMC4) 

CALL REDPI (EKHtS)J 
EKM(6>=DNM(2) 

CALL EKHARD { EJD,DNH, PLA ) 

DO 1 I=l»6 

Y(II=ZERO 

Y<H-6)=DNH(I) 

1 Y(I+12)=ZER0 

YC7)=Y(7)-6NCK£(2 J 
Y(8)=Y(8)-ENCKEf 3) 
y(10>=Y{10)-eNCXE(4) 

HRITE(6,76JY,yp »PAR»tNCKE ,EKH 

EHUM=3012159753997540.0C 
CALLPREITR (00 ,£HUM,EKM, 12,AMM,PM) 

REWIND 1 

CALL OVOPFl (18,Y,YP,K0,EP,RATl,RAT2fHMIN»HMAX,DELT,KST,IHL,PAR 
P,FUNEPH,OUTEPH,KQ,YNN»DT,0> 

WRITEC6,76) ENCKE,AHH,PH 



Program A (Cont) 


CALL SEFODI 11 »2* 18 sTW* 9,3 » FFf PD ) 

REWIND 4 

ET-215.5D0 

OU 33 1=1,53 

DO 32 J=l,15 

READ (2) FD 

DO 31 M=l,18 

DO 31 N=l,3 

31 FH(N,H,J)=FDCN,H) 

32 CDNTI^4UE 
ET=ET+7-DO 
WRITE(4) ET,FH 
BACKSPACE 2 

33 CONTINUE 
REWIND 4 
READt4) ETjFH 
WRITEI6981) ET,FH 
00 34 1=1,50 

34 READ (4) ET 
READ(4) ET,FH 
WRITEt6,81> ET,FH 
STOP 

61 F0RMATJ6F13.6,2X> 

76 F0RHAT(5X,6D19-11) 

81 FDRHAT(//10X,F10.3/,(2X,6D19.1in 
END 
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APPENDIX A/1 


Flowchart of the Simulated Ephemeris Generator 
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APPENDIX A/2 
Sample Input 




MAIN PROGRAM : B 


PROGRAM PURPOSE : To generate optical observations of lunar points. 


INPUT PARAMETERS : Examples of card input are given in 
Appendix B/2. 

A. Simulated lunar ephemeris. 

1. The simulated ephemeris is assumed to be stored 
on a direct access disk device designated unit 4 
(see the main program description for ephemeris gener- 
ation). 


B. 


C. 


Scalars* 

1. STl, ST2 and STS extended precision variables 

are described in the subroutine description OPTOBS. 
Basically, they define the conditions under which a 
point can be observed, 

2. . JOK is a single precision integer variable denoting 

the type of observation (see OPTOBS). 

3. APP is an extended precision variable which contains 
one-half of the field angle in degrees. 

4. ET is the epoch of observation in Julian daj's minus 2440000. 
Vector (extended precision). 

1. X is a 9-element vector input only if the observing 
station is on a satellite which was not the case in the 
example given. For examples of satellite borne observa.’- 
see [Papo, 19711. 


OUTPUT PARAMETERS : 

A. The subroutine computes the simulated optical observations of 

the moon for a bundle of nominally 30 observations of lunar points. 
The angles created are printed (unit 6) and punched (unit 7). 

PROGRAM DESCRIPTION ; A flow chart of the program logic is given in 
Appendix B/1. 


SUBROUTINES REQUIRED : 


A. O. S. U. Project Library: 

1. EPHITL 

2. OPTOBS 
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REFERENCES : 

A. Papo, Haim B. (1971). "Optimal Selenodetic Control, " 

The Ohio State University, Department of Geodetic Science , 
Report No. 156. 

B. Sprague, M. (1961). "An Investigation to Improve Selenodetic 
Control on the Lunar Far Side Using Apollo Mission Trans- 
Earth Photography, " Department of Geodetic Science , 

Report No. 155, June. 
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Program B 


C 


C 


c 


c 


HAIN PROGRAM FOR GENERATING OPTICAL OBSERVATIONS OF THE MOON 


IHPLlClTREAL*aiA-H,0-Zl 

DIMENSION X{9) » BUNDLE <2 »30 ) ,PCSP ( 31 t T (IB ) tTE(3« 3) 
DIMENSION J0KK(2) , I SATT (2 1 . IPOIfU ( 2 1 tUNU ( 2) t AKA (2) 


dimension FH(3.1b»X5J 
COMMON /OPTO/ T8t WTER»SSUN, 
COMMON /EPHEH/ ETE»FH 
REMIND 4 
READ (4> ETE.FH 


ETER,PCSP,ST1,ST2,ST3,IFLAG 


ISET*l 

READ (5 )ST1 1 ST2 f ST3 

1 REA0l5»61fEN0*99) JOK.,APP 
APP - SIN OF FIELD ANGLE / 2 
IF(J0K-91 A,5»99 

4 READ (5,60tEND=99) ET 
CALL EPHITL (ETt4»Y) 

DO 3 M=l»3 
X{H1=Y(H) 

X(H+3)=Y(H*6) 

3 X(M+6)=Y(H+12) 

GOTO 6 

5 R£ADI5,60,EN0=99)ET,X 


6 CALL OPTOSS (£T,X,JOK,APP, BUNDLE) 
If(IFLAG-ll 2,1,2 
2 WRITE(6,70}IS£T,ET 

HRITE(6,75)HTER»SSUN,ETER 

ISET=1SET+1 

IF(J0K.EQ,9) GOTO 7 

R=OSGRT (PCSP(l>»^2fPCSP(2)**2) 

ALL0N=DATAN2 ( PCSP { 2 ) ,PCSP ( 1 n =»57 .295779500 
ALLAT=0aTAN2 (PCSP(3) ,R)*57.2957795D0 


HRI7E(6,72)J0K 
MR1TE(<>,771ALL0N,ALLAT 
GOTO a 

7 WRITE(6,73) 

8 HR1T£(6,76)PCSP 
WRITE(6, 731TB 

WRITE(6,71) (1, (BUNDLE (J, 11,3=1,2) ,1=1,30) 

lSAT=J0X/9 

K0UNT=1 

DO 9 1=1,30 

IFIBUNOLEd,!) -EO.O.OO) GOTO 9 


K0UNT=K0UNT+1 

lH=H0D(KaUNT,2)+l 

IPOINT(IM)=I 

UNUI 1MJ=6UN0LE( I, II 

AKA{IM)=BUNDLE(Z, 1) 

IFnH.NE.2) GOTO 9 

KRITE(7,79) £T, J 0 K,ISAT,IP 0 INT( 1 ) ,UNUm 


,AKA(1),ET,J0K,1SAT,1P0INT 


- 171 - 



Program B fCont) 


P(2)tUNU(2)fAKA{2| 

9 CONTINUE 

IF(1H.EQ,2) GOTO 10 

HRITE{7t79) ET,J0K»1SAT,IP0INT{1) ,UNU{1) *AKA(1) 

10 CONTINUE 

WRITEI7,80) (ET,J0K.lSATt{TB(ItJ)»I=l»3),ISET,J»l«3) 

WRITEl7,80»ET,J0K,lSAT,PCSPfISET 

WRITE(6«S1) X 

GOTO 1 

99 STOP 

60 F0RHAT(5F15.10) 

61 FCRMATt IlOtPlO.5) 

70 FORMATtlHl /10X»»SET NUMBER »t lA^tlOX » ‘EPOCH OF OBSERVATION IN UO 
P - 2A40000.) 1SSF17.9/ ) 

71 FORMAT( 10XfI5,2025.16) 

72 FORHATtlOX, 'STATION OBSERVING IS OBSERVATORY N0',I3/ ) 

73 FQRHATllOX, 'STATION OBSERVING IS ON A SATELLITE'/ > 

74 F0RHAT(2(10X,SD20.10/J) 

75 FORHATtlOX* 'SELENUGRAPHIC LONGITUDE OF TERMINATORS ANO SUBSCLAR PQ 
PINT'/ 10X,3F1C,3/ » 

76 FORMAT ( lOXt ‘SELENOGRAPrflC COORDINATES OF PROJECTION CENTER'/ 10X*3 
PF15.3/ ) 

77 FORHATt 10X» 'TOTAL LIBRATION : IN LONGITUDE* »F8-3, * - IN LATITUDE* 

P,FB.3/) 

78 FORHATf lOXt 'ORTHOGONAL TRANSFORMATION MATRIX FROM THE B1B2B3 TO TH 
PE ECLIPTIC SYSTEM'/3(12X,3D25.16/)/10X*'0PTICAL OBSERVATIONS NU AN 
PO KAPPA'/) 

79 FORHAT{2{F6.2fIlt212,D13.6t016.9) ) 

80 F0RHATtF6.2,Il,I2»3D23.16tI2) 

81 FORHATt/lOXf 'SIMULATED LUNAR EPHEMERIS AND EULERIAN ANGLES * EARTH 
P EULERIAN AN6LFS*/3(12X,3025.16/1) 

END 



APPENDIX B/1 


Flowchart for Generating Optical 
■ Observations of Lunar Points 
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APPENDIX B/2 
Sample Input 


CO 

o 

• 

-0.34 

.34 




3.C07 



3.007 

224.93 



428.87 



2.007 



3.007 

225.01 



430.92 



3.007 




2.007 ■ 

254.92 



431. 



2.007 



3,007' 

255. 



459.93 


** 

3.007 



2.007" 

283.9 



460.01 



2.007“ 


““ 

3.007” 

293.93 



461.98 



3.007 


' 

2.007 

312.88 

2.007 


462.06 

3.007 

312-96 

3.007 


490.96 

2-007 

313.93 

2.007 


491.04 

5.007 

314.01 



519.94 


■ 

3.007 



2.007 

342.92 



520.02 


■ 

2.007 



3.007 

343.0 



548.91 


- . . 

3.007 



2.007 

371.91' 



548.99 



2.007 



‘3.007 

371.99 

3.007 


578.9 

2.007 

400.9 



578.98 
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MAIN PROGRAM : C 


PROGRAM PURPOSE : To adjust the numerically integrated physical librations 
in the simulated environment of CPapo, 1971] to the forced physical 
librations given in [Echardt, 1970], 

INPUT PARAMETERS : An example of card input to the program is given in 
Appendix C/2. 

A. AU through M are input parameters required by the JPL 
subroutines for reading the Development Ephemeris (DE-69) 
on tape as described in [O’ Handley, 1969]. 

B. PARAM is an extended precision vector with 5 elements which 
are described below. 

1, The first element is the initial epoch of integration 
in Julian days minus 2440000. 0 days. 

2, The second element is the final epoch of integration 
' 'in Julian days minus 2440000, 0 days. 

3, The third parameter is the step size (in days) 
used to initialize the numerical integration routine. 

4, The last two parameters are not used as input parameters. 

C. OSH is a 6-parameter vector which is used to input the 
initial values of the physical libration parameters and their 
time derivatives (r^. cr^, p^, r q, p^). 

D. BE is a 6-element vector which is used to input the weight 
matrix of the physical libration parameters. 

E. ALF, BET, GAM, TEQ are inputs of the dimensionless 
ratio of inertia of the moon (a, .8, y) and the inclination 
of the lunar equator to the ecliptic in radians, 
respectively. 

F. The parameters CCE, CCS and P are described in the 
subroutine FUNPL5. 

OUTPUT PARAMETERS 

A. The output of the program is the comparison between the 
simulated solution described in [Papo, 1971] and the actual 
case described in [Eckhardt, 19701, Examples of the output 
are given in [Papo. 1971], 

PROGRAM DESCRIPTION : A flowchart of the program logic is given in AppendixC/l 
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SUBROUTINES REQUIRED : 

A. O. S. U. Project Library: 

1. DVDPF 

2. FUNPL5 

3. OUT AD J 

B. Fortran Scientific Subroutine Package: 

1. DMINV 

2. DGMPRD 


REFERENCES : 

A. Eckhardt, D. H. (1970), '*Lunar Libration Tables, " The Moon, 
VoL 1, No. 2, February. 

B. O’ Handley, Douglas A. et al. (1969). "JPL Development 
Ephemeris Number 69" Technical Report 32-1465. 

C. Papo, Haim B. (1971). "Optimal Selenodetic Control, " 

The Ohio State University, Department of Geodetic Science , 
Report No. 156, 
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Program C 


C MAIN FOR ADJUSTING INTEGRATED NUMERICALLY PHYSICAL LIBRATIONS TO ECKHARDTS 
1MPLICITREAL=»8{A-H,a-Z} 

REAL=w» EP J 60 ) tHHI N, HMAX » EMX ,RAT I ,RAT2 
DIMENSION KD(60) ,KQ(60) ,YN ( 60) ,0T ( 20t60) 

DIMENSION P^3),fiEBC3),SU^^(J),CNuRH(9 ,9 )tCVECT(9' )»G(3t3)t 
•PSEV(6) ,H(13)tH(6,12) ,F(4ViPARAM(3) ,CSH(6) »BE{6) »QI31» 

PY(60),YP(60),LVEC(9 ),MVEC(9 ),TE219 ),C0REL(9 ,9 ),WEC9 ) 
DIHENSIONCNORSCbtB) ,CV£C8 18 J ,LVEtJ(6) ,HVE6 f 6 ) , TEC (3 ) ,CCRE8 ( 8 1 8 } 
COMMON/MAIF/ALF,BETfGAH,TEgfCCt,CCS,P,& /EXE/T00»3EBtSUN ,NUH3 
P/MAICUT/CNORH,CVECT /ALL/ SEV .ROM ,NUT, KLU 

P/CETBLl/AUtREA,TPDt£MR/C£TBL2/ICH,ICE,H/CETELA/rifF 
DATA ONEfTWOfFOR/l,DO,2.COtA,DO/ 

EXTERNAL FUNPLS .QUTAOJ 

REA0(5,60)AU»REAtTPD»EKR 

REA0(5t65)1CE,M 

REA0{5,60)PARAM 

READ (5 ,65 >>40 

READ<5,60)HE 

READtS, 60) HMIN, HMAX, DELT, RATI, RAT2 

REWIND I 

1CH=1 

NUT=2 

KLU=1 

2 CONTINUE 
REAO(5,60)0SH 
READ{5,60)Q,TE0 
READ(5,62)CCE,CCS 
REA0(5,62)P 

• ICHE=1 

101 ICHE=1CHE+1 

1 Al=ONE/(Ot31-(TWO-FOR*Q(2))=»0(l)) 

81=0NE/( (0NE+QI2) )*Q (3 )-{ TWQ-TH0*Q{2 1 )*0C1 ) ) 

ALF=QU)*{Q{3)+TWC»0t 1))*A1 
BET=Q(2) 

GAM=-FOR*5{2)*Q(l )*B1 

MR1TE(6,79) PARAH,HE ,Q , ALF ,BET,GAH,TEq,CC£,CCS,P 

ROMsO.OO 

NUMB=0 

16 DO 3 1=1,6 

3 YCI)=OSH(I) 

IF(KLU.GT.l) GOTO 8 

17 A1=AI»A1 
B1=B1*B1 

G(l,l)= F0R^QI2)*(0NE-Q{2) )*Q(3)*A1 

G(l,2) = {a(3)**2-F0R*Q(1)=^+2>»A1 

GU,3)=-F0R*D(2)»(0NE-0(2))=ytl)*Al 

G(2,I)=0.00 

G(2,2)=1.D0 

G(2,3)=0.00 

G(3,1)=-F0R*LM2)*<0NE+Q(2))*Q(3)*B1 
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Program C (Cont) 


G{3,2)=-F0R*0C 1>*C0C3)-TW0*QC 1) )^E1 
G(3»BJ:= FOR*Q{2)^{QUB^Qi2))^Oil)^Sl 
C PARTIALS OF ( ALF, BET, GAM ) VS* (C22, BETt C20 ) 

DO 9 K=l,9 
CVECTCK)=O.DO 
DO 9 1=1,9 
9 CN0RH(K»n=0*D0 

18 DO 5 1=7, NQ 

5 YCI)=0*DO 

00 6 I=7,42t7 

6 Y{I)=1.D0 

C INITIALIZATION OF TRANSITION AND SENSITIVITY MAIRICES 

8.T00=PARAHJ D + l.ODO 

CALL DVDPF (NQ,Y , YP,KD , EP , RATI ,RAT2 , HHIN,HHAX , CELT, KSTtr 
PIHL,PARAH,FUNPL5,0UTADJ,KQ,YN,DT) 

HRITE(6,79) ROM 
IF( 1CHE*EQ*3) GOTO 998 
DO 14 1=1,8 
CVEC8(I)=CV£CT(I) 

DO 14 J=l,8 

14 CN0R8CI,J)=CN0RH( I,J) * 

CALL DHINV (CNGR8 ,8 ,OET, LVE8 ,MVE8 ) 

CALL DGHPRO { CN0R8,CVEC8 , IE 8, B , 8 , 1) 

DO 10 1=1,6 

10 0SHCI}=0SH(I)-TE6(I) 

WRITE(6,73) TE8,CNOR8 
Q{l)=Q{l)-TeS(7) 

Qt2)=Q(2)-TE8{3) 

DO 11 1=1,8 
DO 11 J=l,8 

II- C0RE8{'I,J)=CN0R8(I,J)/ OSORTt CN0R8{ I , I ) ♦CNORS I J, J ) ) 

HRXTE<6,78) CORES 
IFCICHE*EQ*1) GOTO 101 
NQ=6 

OELT=10*D0 
KLU=2 
GOTO 101 

19 N0=60 
KLU=l 
GOTO 101 

998 STOP 

60 F0RMATC6F13*5) 

62 F0RHAT(4D20*13> 

65 F0RHATU6I5) 

78 F0RHAT(/<5X,S013.5n 

79 FORMAT (5Xt6DIB.10) 

80 F0RHAT(5X,9 011*4) 

81 FORMATtSX, 8011-4) 

90 F0RHAT(5D16*9) 

END 
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APPENDIX C/1 
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APPENDIX C/2 
Sample Input 


— K9597893*; 6378 ^1492 * !• 


81.301 


570. 


1 . 


— 203 r 

-.00001 A61873. 023799910601.00001 


— 1000 / 
.0003968 
—.8926626 


10 . 

.0006268 
0 16 


1000 . 

.00023 


“ ' O.OOCOl ' 

8879268.000008645041 

' ■ 1000 . 

.026769 


1 1 

-.00021278188 

loo; ■■ ' ' 


-.00016866572 

1000 . 
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MAIN PROGRAM ; D 


PROGRAM PURPOSE : To process earth-based optical observations of elected 
lunar points in a weighted least squares procedure. 

INPUT PARAMETERS : An example of the card input required for the program 
is given in Appendix D/2 (except for estimates of the TQ). 

A. PARAM is an extended precision vector of 5 elements. 

1. The first element is the initial epoch in Julian days 
minus 2440000. 0 days. 

2. The second element is the final epoch in Julian days 
minus 2440000. 0 days. 

3. The third element is the initial step size in days. 

4. The remaining elements are not used for input. 

B. NQ is a single precision integer used to input the number of 
integrands. 

C. OSH is an extended precision 6 element vector used to input 

the initial values of the physical librations (r, <j, p) in radians 'sr. 
their time rates (r , a, p) in radians per day, 

D. The vector SEV is described in the subroutine FUNPL5 descript:o- 

E. The variables HMIN, HMAX, BELT, RATI, RAT2 are describt:d 
• in the subroutine DVDPF2 description. 

F. The variables Q, TEQ, CCE, CCS, P and G are described 
in the subroutine FUNPL5 description. 

G. TJ is an extended precision array of 3 x 22 elements which 
contain approximate (initial) selenocentric positions of the 
lunar points in kilometers. 

H. JOBN is a single precision integer vector of 10 elements 
which is used to input the identification number of batches 
of optical observations to be processed. 

I. Simulated ephemeris data is assumed to be available on a direct 
access disk device (designated unit 4), The parameters read 
from the disk are ETE and the array FH (see EPHITL subroutine 
description). 

OUTPUT PARAMETERS ; 

A. Basically output is created in the subroutine OUTADJ which 
in turn is called by DVDPF2 Samples of the 
normal equation are given in f.Papo, 1971]. 

PROGRAM DESCRIPTION : The program was formulated according to the theoreti' ■ 
development given in fPapo, 1971] Chapter 2. A basic flowchart 
is given in Appendix D/1. 
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SUBROUTINES REQUIHED : 

A. O, S. U. Project Library: 

1. DVDPF2 

2. FUNPL6 

3. OUTGAJ 

REFERENCES : 

A. Papo, Haim B. (1971). "Optical Selenodetic Control, " 

The Ohio State University, Department of Geodetic Science , 
Report No, 156. 
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Program D 


C HAIM FOR ADJUSTING SIMULATED EARTH-BOUND OPTICAL OBSERVATIONS 
IHPLICITREAL^a (A-H,0-2 J 
REAL^4 EP(60) ,HM.IN,HHAX,£HX,RAT1,RAT2 

DIMENSION K0C60 J tKQCSC) f YN(60) ,0T(20,60) »FH<3tX8tl5>,VEC W5) 
OIMENSIONPOSC3) fEH0T<3,3) ,ENU<22) t£KA(22) ,TJ(3*22) »SIG6(3,3> 
P»SIGZ(3»3)flP0(22) t JOBNdO) 

DIMENSI0NP(3),(iED(3) ,SUNC3) iG(3t3NY(60) ,YPf60)» 

PSeV(6) ,H(13)tH<6»12) ,F{A)tPARAM(5 ),CSH(6) tQ(3) 
COHHON/MAIF/ALFTBETtGAMtTEQ,CCE,CCStPtG /EXE/TCO t SEB tSUN ♦NUMB 
P/HAICrjT/PO$*£HBT,ENU,EKAfTJ,SIG6^SIG7,SIG2»K3»K2»IPO,IST,IlCW 
P /ALL/SEV,ROM,NJT,KLUTN2tKMf ISKiP /EPMEH/ETE»FH 

DATA ONEtTWO»FCR/l*DOf2*DO,4*UO/ 

■ DATA PI /3*141S92fe53589793D0/ 

EXTERNAL FUNPL6 rOUTGA J 

READ{5»60)PARAM 

REAO(5,65)NQ 

READ{5»60JOSH 

REA0(5»603SEV 

REA0(5T60)HMINtHMAX,0ELTTRATltRAT2 
READ(5f 60)Q»TE0 
READ{5f62)CCEtCCS 
REAOt5f62)P 
REA0(5f61) TJ 
READ 15, 65) JOBN 
KLU==1 
NUHB=l 
KH=30 
1SK1P=4 
I1CH=1 
.111 = 1 
DO 4 . 1 = 1,3 
DO 2 J=lt3 
SIG6(I ,J)=O^DO 
2 S1G2(I ,J)=0.00 
SIG6(I,1)=.0200 
4 SIG2(I,I)=*0400 
SIGZ(1,1)=*004DO 
SIG7=-0i00 
REWIND 8 

709 JOB = JOBN nil) 

1F(J03*EQ.99) GOTO 998 

REWIND 3 

REWIND 4 

READ (4) ETE,FH 

NUT=0 

N2=44 

K3=l 

K2=2 

REWIND K3 
00 28 1 = 1,151 
00 21 J=l,75 
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Program D (Cont) 


C 


c 


21 VEC(J)-0.00 
WRITE(K3) VEC 
28 CONTINUE 

1 A1^0NE/CQ(31-(TW0-F0R*Q(2) )^QC1) ) 

Bl=0NE/t(CNE+Q<2) )*C«3)-tTWO-TWO*0(2) J*Qm) 

ALF=Q{2J*(Q<3)+TWQ«'Qa))*Al 

BET=Qt2) 

GAK=-f0R*0{2)»Q{U>i'Bl 

WRITE{6»79> SEVt PARAH ,Q,ALF,BETtGAM-tTEQtCCE»CCStP 

16 DO 3 1=1»6 
3 Y(I)=OSH(Il 

17 A1*A1*A1 
B1=B1*B1 

GUtD- FaR*Q{2)*<0NE-Q(2) )»0{3)*A1 

G(l»2)={Q(3J*»2-F0R*0tl)*^2)»Al 

GC1.3J=-FOR*OI2I=*(ONEHS(2))»QUJ*A1 

Gt2,l)=0.00 

G(2,2)=l.DO 

GC2i3)=0,U0 

6(3,l)=-FOR*Q{2 J*{0NE+Q{2) )*Q(3)*81 
GI3i2)=-P0R*0(l )* (0(3 )-THG*Q( 1) )*B1 
G(3t3J= FOR*Q(2)*(ONE+Q(2) J *0(11*01 
PARTIALS OF < ALFt6ET,GAH> VS. (C22«BET»C20) 


10 DO 5 l=7tNQ 

5 Y(1)=0,D0 

DO 6 I=7,A-2»7 

6 Ytl)=1.00 
15 CONTINUE 

INITIALIZATION OF TRANSITION AND SENSITIVITY MATRICES 
8 TOO=PARAH(1)+1.0DO 
DELT=1.00 

CALL DVDPF2(N0»Y ,YP tKO»EP t RAT 1 1 RAT2 tHMINt'HHAX»OELTtKSTf 
PIHL, PARAH, FUNPL6.0UTGAJ ,KQ,YN,DT,1) 

REWIND K3 
DO 151 1=1,151 
REA0(K3) VEC 
151 WRITE(B) JOB, VEC 

iii=m+i 

DO 801 1=1,3 
801 SIG6I1,! J=.2D0 
S1G7=.100 
GOTO 709 
998 STOP 

60 FORMAT (6F13. 5) 

61 FCRHAT(2{10X,3F10,A)J 

62 F0RHAT(AO20.13) 

65 format (161 5) 

79 F0RHAT(/(5X, 6018. ion. 

END 
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Flowchart of the Optical Data Processor 
(Earth-based Observations Assumed) 




INITIALIZED INTEGERS* 

KLU- Indicates if portiol 

dervJotive ore integrate 
or read from disk 

NUMB - Not used 

*KM- Number of optical 

bundles to be processe 

ISKIP- OUTGAJ 

IICN - Bundle sequence 
number 

111- Job sEriai number 

SET INTERGERS 

nut- If equol to 2 , lenntnQ 

NZ - Two times the number 
lunar points 

K2- K3 - Unit numbers of 
auxiliary device 


135 













APPENDIX D/2 


Sample Input 


222.5 579. .5 .000001 

• 60 ' ‘ 

-.0000607 18 23-.00&81857^A5.0C06061C22«6. 000027946591. 00530^180A4 .000062951555 

1.196721511 0.08051275 .2299715022 -.00092921392 

.005 .5 1. .000001 .00001 

.0000207 '.000629 -.000207 .026769 

.892668 D 16 

31 32 99 
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MAIN PROGRAM : E 


PROGRAM PURPOSE : To create simulated laser ranges. 

INPUT PARAMETERS : 

A. Integers (single precision) card input. 

1. IDP is an 8-element vector used to input earth 
= station identification numbers. 

IDL is an 8 -element vector used to input lunar 
station identification numbers. 

B. Matrices (extended precision) card input. 

1. XP is an 8 X3 array containing the U, V, W geocentric 
coordinates of the earth stations in kilometers. 

2. XL is an 8 X3 array containing the X, Y, Z seieno- 
centric coordinates of the lunar reflectors in kilometers. 

C. Scalar (extended precision) disk input, 

1. ETE is the Julian day epoch for the lunar ephemeris 
stored on disk. The parameter is described in the 
EPHITL subroutine description. 

D. Matrix. (extended precision) disk input, 

1. FH is a 3 Xl8 Xl5 array containing tabulated ephemeris 
quantities as described in the EPHITL subroutine 
description, ■ 

OUTPUT PARAMETERS : The program prints and punches output for the 
simulated laser ranges. The parameters are: 

A. Integers (single precision). 

1. IDl is the station identification number of the observing 
station. 

2. ID2 is the reflector identification number observed. 

B. Scalars (extended precision). 

1. JED is the Julian date of observation minus 240000. 0 
days. 

2. D is the simulated observed distance in kilometers, 

3. ZD is the zenith distance of the observation in degrees, 

C. Vector (extended precision). 

1. X is the topocentric Cartesian coordinates of the 
reflector as described in the LADIS subroutine 
documentation. 
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PROGRAM DESCRIPTION : A flowchart is given in Appendix E/l to outline 

the programming of the theory developed in [Fajemirokun, 1971]. 
Sample input is illustrated in Appendix E/2. 


SUBROUTINES REQUIRED : 

O.S.U. Project Library: 
1. EPHITL 
. 2. PMAT 
3. LADIS 


REFERENCES : 

Fajemirokun, F. A. (1971). "Application of New Observational 
^sterns for Selenodetic Control, " The Ohio State University, 
Department of Geodetic Science. Report No. 157. 
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Program E ’ 


C 


implicit realms (A-H»0-Z) 

MAIM PROGRAM PGR COMPUTING SIMULATED LASER DISTANCES 


DOUBLE PRECISION JED 

01MENSi0NXP(8,3> tXU6t31»XEl3J »XM(3Ti XC5I3»,P£( 
*FH{3»iSf 15) rUU( 16) ♦lOH(fc) ,1DLIU 
COMHON/EPHEM/ETEtPh 


3),PHI3»3) »XC3)» 


I*l 

READ{5.8inOPU).tXPU.J)»J=L»3» 

WRIT El 6,81) IDPI n»(XPII,J),J=l,3) 

30 CONTINUE 

00401=1,3 

READ|5,8l)IDLll) , IXLf I ,J) , J=1 »3) 

WRITE C 6,81 ) lOLI 1 > fiXwl i»J)'f J=l,3) 

40 COKTINUfc 
REWIND 4 
READ! 4) ETh,FH 
JtD=222.5D0 

45 CALLEPHiTL(J6D,4,wU) 

D050I=1,3 

xcetu=tuu) 

50 CONTINUE 
PHlM=Cum 
PS IK=0UI6) 

THB1AM=QUIV) 

CALLPHATIPSIM,THETAM,PHiH,PM) 

PHIE=QU(13> 

PSIE=0U(14) 

THP1AF=QU115) 

CALLPHATIPSl£,TH£TAt,PHiE,P£) 

00601=1,5 

102=IDL(1) 

XMIl)=XLU,i) 

XM(2)=XL{ l.,2) 

XMI3)=XLI 1,3) 

J=1 

101=IDP(J) 

XbU>=XPlJ,i) 

XE(2)=XPtJ,2) 

Xei3)=XPI J,3) 

CALLLAOISI JE0,P£,PH,XE,XM,XC£,X,0,ZD) 

lPC0.£C.O,0DO)l>0T060 
WRITE lt>,91) JED, l!ji,lD2,Xt[>,Z0 
WRITE (7, 92 )JE0,i0l,ID2,X,3,20 
60 CONTINUE 

JE0=JeD*0.5D0 

lFIJED.G£.250.000)GUr070 

G0T045 

70 HR1T6|6,93) 

STOP 
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Program E (Cont) 


oi FCHMaTI i5j3D^0«9) 

l-URMAr(i2,t>Xf JFiO«‘^) 

VI |-(jRMAT(/t2X,F7-2,2i5» 3019.10) 

V2 F0RMAnF7-3,212,Vhio.7,f3.i) 

V3 FORMATC »1« tiOXv *JOb Rt^^UlRtD IS DO^ifc*) 
Vo FORMAT {3D25. 12) 
tND 
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APPENDIX E/1 


Flowchart for Laser Range Generation 
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iM m 


APPENDIX E/2 
Sample Input 


0 . 1-33 

1 0. 1391800000+04 

XT. 1 6 5 ZCKTOrrcnr+TTzr 
0. 1554000000+04 


0 A01 ;<7 ?0000+04 0 . 1. 93 560OQ0D+ 0 2 APOLLO 11 

=D'75^0 000001)+ ( 0^ -0. llUOOOOOOlJ+lB APULLl' -1^ 

'">.9920000000+02 0.7620000000+03 APOLLO 15 
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MAIN PROGRAM : F 


PROGRAM PURPOSE : To adjust simulated laser distances to the 
moon. 

INPUT PARAMETERS : 

A. Integers (single precision). See Appendix F/2, 

1. NE is the number of earth station coordinates. 

2. NM is the number of lunar reflector coordinates. 

3. NV is the number of unknowns in the solution. 

4. NOB is the number of observations in the solution. 

OUTPUT PARAMETERS : Output for the program is accomplished in the 
subroutine LASOLV. 

PROGRAM DESCRIPTION : A flow chart of the program is given in 
Appendix F/i. 

SUBROUTINES REQUIRED ; 

O. S. U. Project Library : 

1, LASOLV 

REFERENCES : 

Fajemirokun, F. A. (1971). "Applications of New Observational 
Systems for Selenodetic Control, " The Ohio State University, 
Department of Geodetic Science, Report No. 157. 
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o u 


Program F 


MAIN (1ST) PROGRAM FQrt AbJUSTiNG SlMULATtD LASER nlSTANCES 
i/.PLlCl IREAL^b U-H ,0-Z ) 

DiMENSiUN lDP(B),XP{b,3)»iDL(t),XLtiJ*b)TB(2iO»63)7EMI\V(2iO)f 
^wLA(63,63) ,PX(o3) tT£MP(2iO) ,LVcC(63) *HVEC{fe3) ,EN(63 ,o3) > 

=tCuRtL( 63t 63 ) 

bIHENSIQN bM(63,2iO) ,W (210) 

1 CONTINUE 
DU15I=1»210 
D015J=1*63 
8a,J)=0.000 
b.M(J, i> =0.000 
13 CONTINUE 

DGZOI=lf2iO 

EHINV( 1 )=1. DC/ (0.1300=^0. 1300) 

WU) =0.000 
20 CONTINUE 

KtADi 5,10)NErNM,NUtN0b 
WR 1 TE ( 6 , 10 ) NE tNM t NU, Nl tj 
NP=NE/3 
NL=NH/3 

call LASOLVl IDP7XP»iDL,XL»b,LMiiJV,QLA»PX»TEMP»LVECtWVEC,EN» 

=?LCRE‘L»NE»NH,Nb7NDb»NPtI'iLTi}H*ri) 

STOP 

10 FORMAT (A 15) 

LhD 
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APPENDIX F/1 
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APPENDIX F/2 
Sample Input 
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